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ABSTRACT 



The ability of a submarine to maintain ordered depth, especially during periscope 
depth operations at low speeds, is vital for the vessel to perform its mission and avoid 
detection. Modem submarines exhibit an inherent phenomenon that produces an 
undesirable ship response at low speeds, commonly referred to as dive plane reversal. The 
physical parameters that govern this occurrence are related in this thesis to the problem of 
multiple steady state solutions in the vertical plane. 

Generic solution branching, in the form of pitchfork bifurcations, can occur when the 
nominal level flight path loses its stability. A systematic study reveals the existence of a 
critical Froude number, based on the vessel's speed and metacentric height, where this 
branching occurs. Bifurcation theory techniques and numerical computations are utilized 
to classify the effect that geometric parameters, trim and ballast conditions, and 
hydrodynamic properties have on the existence of these multiple solutions 
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I. INTRODUCTION 



One of the most critical functions of a submersible is accurate depth keeping at the 
commanded depth. Such a function can be carried out either manually or automatically, 
especially in cases where human intervention is not possible or desirable Due to the 
technological significance and numerous scientific applications of submersible vehicle 
systems, design issues of appropriate depth keeping control laws have received wide 
attention in the past Such control system designs include linear and nonlinear controllers 
[Refs 1&2], model based compensators [Ref 3], adaptive control [Ref. 4], and sliding 
mode control laws [Refs. 5&6], 

Response accuracy and stability are primary considerations in designing depth 
keeping control law Of paramount importance, in this area, are the robustness properties 
of the particular design; i.e., its ability to maintain accuracy and stability in the presence of 
incomplete sensor and environmental information, as well as actual/mathematical model 
mismatch. The scope of the work in this thesis is to demonstrate a potential loss of 
stability that may occur when a submarine is operating at low speeds. This loss of stability 
can occur regardless of the particular means used for depth control. The study is 
accomplished through the use of an eigenvalue analysis, a steady state analysis and a 
controllability analysis [Ref. 7], It is shown that such a loss of stability is accompanied by a 
slow divergence of trajectories away from the commanded path Solution branching 
occurs in the form of generic pitchfork bifurcations [Refs. 8, 9], A complete 
characterization of the problem is given utilizing singularity techniques, which have been 
proven to be very useful in the analysis of similar problems [Refs. 10, 11, 12], The use of 



bifurcation theory allows the crucial vehicle parameters that govern the problem of 
solution branching to be determined, and the develop guidelines to prevent its occurrence. 

Finally, a new look at the problem of dive plane reversal [Ref. 13], based on solution 
branching results is presented The term dive plane reversal refers to a well-known 
phenomenon in submarine operations where, during low speed depth keeping, there is a 
need to reverse the direction of dive plane deflection in order to execute a given change in 
depth Physically, this can be explained by considering the relative magnitude of the 
hydrodynamic forces At moderate and high speeds, the normal force on the submarine's 
hull due to the angle of attack exceeds the normal dive plane force and the boat responds 
to ordered dive plane angles as expected. The phenomenon of dive plane reversal occurs 
at speeds below a certain critical speed in which the normal hull force is less than the 
normal dive plane force and the response of the boat is reversed. [Ref. 14] Vehicle 
modeling in this work follows standard notation [Ref. 15] and numerical results are 
presented for the DARPA SUBOFF model [Ref. 16] for which a set of hydrodynamic 
coefficients and geometric properties is available. Special emphasis is given to identifying 
the proper non-dimensional parameters in the problem, so that extension of these results 
to full scale models and other designs is possible using minimal experimental and/or 
analytical results. 
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II. VEHICLE MODELING 



A. EQUATIONS OF MOTION 

1. Introduction 

For the purpose of this problem, and the subject of the maneuvering and motion 
control of the vehicle, the following assumptions are made: 

1 ) The vehicle behaves as a rigid body; 

2) The earth's rotation is negligible as far as acceleration components of the 
center of mass are concerned; 

3) The primary forces that act on the vehicle have inertial, gravitational, 
hydrostatic and hydrodynamic origins. 

2. Coordinate Systems And Positional Definitions 

A global navigation frame, OXYZ, bs shown in Figure 2 1, is defined with origin, 
O, and a set of axes aligned with directions North, East and Down. This produces a right- 
hand reference frame with unit vectors 7, J, and K Ignoring the earth's rotation rate in 
comparison to the angular rates produced by the vessel's motion, it can be said that the I , 
J, and K coordinate frame is an inertial reference frame in which Newton's Laws of 
Motion will be valid. As seen in Figure 2 1, a vehicle's position R g . , in this frame will have 
the vector components, R . = [ XJ + Y 0 J + Z 0 K ] . A standard convention that is used 

in both aircraft and marine vehicle dynamics places the Y axis to the right while looking 
along the X axis, and the Z axis is positive downwards. Figure 2.1 also shows a vehicle 
with some general attitude in the original coordinate system. 
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Figure 2-1. Coordinate Axes Convention. 
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Now define a body fixed frame of reference Om, with the origin O , located on 
the vehicle centerline, moving and rotating with the vehicle, in which the vehicle's center 
of mass, G, has some position other than the origin of the vehicle fixed frame (refer to 
Figure 2-1). This origin will be the point about which all vehicle body force will be 
computed in later sections of this chapter. The convention used for the vehicle fixed 
frame, with unit vectors T , ] , k , is that its origin lies at the ship's center (origin at the half 
length), the x axis at main deck level, parallel to the longitudinal centerline, with the r axis 
vertically down The vessel's center of gravity and center of buoyancy do not generally lie 
at the origin of the body fixed frame, nor are they collocated. These points are denoted by 
G and B respectively. The vector components of p G and p B are thus [ x G J + y G ] + z G k ] , 
and f xj + y B ] + : B k J The location of the center of mass is important because 
Newton's Laws of Motion equate forces on a body to the rate of change of linear 
momentum of the center of mass and moments about the body center of mass to the rate 
of change of angular momentum. This is a particularly important point to bear in mind for 
ship and submarine motion dynamics as the center of buoyancy is determined by the shape 
of the submerged portion of the ship body while the center of gravity is determined by the 
distribution of the weight over the entire ship body. Having defined a coordinate system 
that will be used to describe a ship's position, there is a need to define angular orientation 
of the ship's body, leading to the definition of translational and rotational velocities and 
accelerations. 

3, Angular Position In The Global Reference Frame 

There are several different ways that the attitude of a vehicle can be described in 
reference to the global frame. The most usual and common method is to define three 
angles, called Euler angles that uniquely define the angular orientation of the vehicle 
reference frame, relative to the global reference frame One problem with the Euler angle 
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approach is that a singularity exists when one of the angles reaches 90 degrees (an aircraft 
in pure vertical flight cannot distinguish its azimuth angle from its roll angle). This 
limitation which can sometimes, although rarely, cause trouble in flight simulations and 
control computations can be overcome by the use of quaternions which introduce four 
rather than three variables to describe angular position In this presentation, however, we 
will use Euler angles as it is the most widely used method, although the use of quaternions 
has found favor in robotics applications and computer graphics. 

While any consistent definition of three base angles would be sufficient, the most 
convenient formulation of these angles is in the widely used military terms of an azimuth, 
elevation, and spin notation The rates of change of these angles do not generally 
correspond to the other commonly used angular rates describing angular velocity of a 
vehicle, that is yaw rate, pitch rate and roll rate except, as will be seen, where motion is 
limited to small angle rotations. 

4. Rotational Transformations 

Vehicle attitude is important in defining position. Under dynamic conditions, the 
pitch or roll can cause problems for manned vessels and the heading is critical in 
navigation. For the purpose of considering angular rotations, consider a forward 
transformation from a coordinate triad aligned with the global reference frame and perform 
a sequence of three rotations to finally align the result with a frame that is assumed to be 
parallel to the vehicle body coordinate axes, and moving with the vehicle at all times. 
Begin by defining an azimuth rotation, \\i , as a positive rotation about the global Z axis. 

Next define a subsequent rotation 0 , (positive up) about the new Y axis, followed by a 
positive rotation 0 , about the new X axis. The triple rotational transformation in terms of 

these three angles, is then sufficient to describe the angular orientation of the vehicle at any 
time. It follows that any position vector, R 0 , in an original reference frame given by 
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R o = fx 0 ,y 0 ,z o J , will have different coordinates in a rotated frame when an azimuth 
rotation by angle y , is made about the global Z axis. 

If the new position is defined by, = [x v y v z x j, it can be seen that there is a 
relation between the vector's coordinates in the new reference frame with those that it had 
in the old reference frame. Using trigonometrical relationships and Figure 2.2 as a guide, it 
follows that, 

X \ = X a cosy + Y a sin y (2.1) 

= -X 0 sin\\t + Y 0 cos\y . (2.2) 

This relation can be expressed in matrix form by the rotation matrix operation, 

< 2J > 

where the rotation matrix \T y2 \ represents an orthogonal transformation. 
Premultiplication of this rotation matrix with any vector, R a , will produce the components 
of the same vector in the rotated coordinate frame. Continuing with the series of rotations 
results in a combined total rotation transformation, 

T(W,v)=T(VT(Q)T(y). (2.4) 

In expanded notation equation 2.4 takes the form; 

cos y cos 0 sin y cos 0 - sin 0 

cosy sinQsinty-sin\\f cos <}> s7«ys/>;0sfr/<() + cosy cos(|) cosQsinty , 
cosy s/z/0 cos <}) + s/>?y s/>7 <|> s/// y s/>?0 cos (})- cosy sin§ cos 0 cos (p 

and it can be said that any position vector in an original reference frame may be 

expressed in a rotated frame with coordinates given by the operation, 

JL = [7m>trM*. (2 5) 
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Figure 2-2. Azimuth Rotation. 
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5. Kinematics 



Kinematics deals with the relationships of motion quantities regardless of the 
forces induced by their prescribed motions. Descriptions of a ship's position both 
translational and rotational, will need to be related to velocities, both translational and 
rotational, prior to extending the situation to accelerations. The connection between 
translational velocity and the rate of change of translational position is straight forward. 
Define the global velocity as. 



R = 



X 

Y 

Z 



( 2 . 6 ) 



This vector will have components that are different when seen in a body fixed frame 
Define the body fixed components of the global velocity vector as [u, v, w]'. These 
components, in terms of the above global quantities are given by the forward 
transformation defined earlier to be, 



u 




~x 


V 

Vi’ 




Y 

Z 



(2.7) 



It is a simple reverse coordinate transformation from body fixed to global coordinates to 
see that, 



X 




U 


Y 

Z 




V 

Vi’ 



This shows that the progression of a vehicle in a global frame clearly depends on its local 
velocity components and its attitude. Put simply, u corresponding to a vehicle's forward 
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speed (surge), v corresponding to a side slip velocity (sway), and w corresponding to any 
velocity component in the local Z direction (heave) and the vehicle's global velocity 
components depend on heading, pitch, and roll attitude. 

The connection between angular attitude and angular velocity is not quite so 
apparent At first sight, it is tempting to define the instantaneous angular velocity of the 
vehicle simply as the rate of change of its angular position defined by the Euler angles. 
This is erroneous however, because the rotation 0 , was defined as a rotation about the 
intermediate frame after a rotation \]i had been made. Vehicle inertial angular rates are 
defined in terms of components that have angular velocities about the global axes. It is 
necessary to relate both Euler angles and their rates of change to angular velocity 
components about the global axes to their components lying along the body fixed axes in 
any attitude The prime reason for this is that it is difficult to construct physical sensors to 
measure rates of change of Euler angles. However, rate gyros in common use today are 
easily constructed to measure the components of the inertial angular velocity of a vehicle 
that lie along the vehicle's body axes. It follows that the instantaneous angular velocity of 
the vehicle can be related to the instantaneous rate of change of angular orientation only 
after considerations of the intermediate transformations used. In other words, if a is 
defined as the angular attitude vector, a = [\|f,0,<t>] ; and the inertial components of the 
vehicle angular rate lying along the body axes as to = [p, q,r ] ; then, a - f( (a). The details 

of the nonlinear functional relations involved are provided by viewing the rate of change of 
the rotation y as a vector quantity lying along the original Z axis. The rate of change of 
the angle 0 is viewed as a vector quantity lying along the Y axis of the first intermediate 
frame, and the rate of change of the angle <)) is viewed as a vector lying along the X axis of 
the final (body fixed) frame. Each of these component rates of change of angular position 
has component parts that project onto the final frame and it is the sum total of all the 
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components that give the total angular velocity as seen in the final frame of reference 
Using the required transformations for the rate components from each Euler angle, we get. 



V 




'o' 




'o' 




V 


<7 

r 




0 

y. 


+ m)T(Q) 


0 

_0_ 


+m) 


0 

0 



with the result; 



( 2 . 9 ) 



p 




-vj/5/>t0 + <j) 


<7 


~ 


vj/ sin 0 + 0 cos (j) 


r 




\j/ co50 cos<{) - 0 sin <{> 



( 2 . 10 ) 



Inverting equation 2.10 yields a solution for the rates of change of the Euler angles in 
terms of the body fixed components of the angular velocity vector, 



V 




p + qsin§tanQ + rcos§tanQ 


0 


- 


qcos§-r sin§ 


y_ 




(qsin <{> + r co 5 <() ) / CO 50 



Notice that for small angular rotations, as expected. 



<\> = p; 

Q = q; 

=r. 

At this point the kinematic relationships between velocities, as seen in the body fixed 
frame, and the rates of change of global positions and Euler angles have been defined. The 
resulting set of differential equations forms a consistent set in that given a set of vehicle 
velocity data versus time, its position and attitude may be computed 
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6. Translational Equations Of Motion 

The global acceleration of the center of mass is derived by differentiating the 
velocity vector, taking into account that the center of mass lies in a rotating reference 
frame Considering the total differentiation, the global acceleration of the center of mass 
becomes, 

Rc = v + djxp G +d>xcoxp G + cbx v, (2.12) 

where v = R o . The translational equation of motion is found by equating this acceleration 
to the net sum of all forces acting on the vehicle in three degrees of freedom (X,Y,Z) One 
important factor to recognize is that the equation of motion derived in this manner is a 
vector equation with the components expressed in the body fixed frame and unit vectors 
T, j and k . This has been deliberately done because the dominant forces acting on a 
submerged body in motion are developed in relation to the shape of the vehicle and are 
more conveniently expressed in relation to the body axes. Equating the applied force 
vector to the acceleration, results in, 

F = m{v + wxp G + d)Xwxp G + (bxp G y . (2.13) 

The applied force vector is composed of gravitational (weight) and hydrostatic (buoyancy) 
forces together with any hydrodynamic forces arising from relative motion between the 
vehicle and the ocean water particles. These will be discussed in further detail in following 
sections 

7. Rotational Equations Of Motion 

To develop the rotational equations of motion, the sum of applied moments 
about the vehicle’s center of mass is equated to the rate of change of angular momentum of 
the vehicle about its center of mass. This equating will provide the necessary equations of 
motion for the remaining three degrees of freedom. In the practical case of marine 
vehicles, however, the statement just made is modified slightly because it is much more 
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difficult to assess the vehicle's mass moments of inertia about its center of gravity (CG) 
As the CG changes with loading; it becomes simpler to evaluate the mass moments of 
inertia about the body fixed frame which lies along the axes of symmetry of the vehicle in 
most cases. The mass moments of inertia for the vehicle to be computed are; 




/ 


1 


xy 


X Z 


I 


I 


>v 


y: 




I :: 



(2.14) 



The angular momentum of the body is thus; 

H 0 = I 0 «>, (2.15) 

resulting in the total applied moments about the origin given by, 

= H 0 +Pg x (2 16) 

Since the rate of change of angular momentum is given by, 

R o = I o G) + G>xH 0 , (2.17) 

and the acceleration of the global position vector is given by, 

R o = v + (bxv, (2 18) 

then the rotational equation of motion in vector terms is given by, 

M a = / 0 dj + (bx( / o d ) ) + m{p G x £ + p G xcbx v}. (2.19) 

The applied moments about the body fixed frame arise from a static balance of 
weight and buoyancy effects to achieve the proper trim and heel, and hydrodynamic 
moments from the forces applied through appendages such as control fins, hydrodynamic 
effects from waves, and hydrodynamic effects from relative motion between the vehicle 
and the water. 

At this point, there are three translational equations obtained from equation 
2.13, three rotational equations obtained from equation 2.19, and six unknown velocities 
(co and v). In itself this is a consistent set of dynamic equations if the weight and buoyancy 
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terms are always self canceling. However, with weight and buoyancy being applied in a 
global vertical direction, and also being applied at different locations, they represent forces 
that are dependent on the attitude of the vehicle. Notice that without weight and buoyancy 
forces, the six consistent equations expressed in the body fixed frame of the vehicle could 
be solved for the vehicle's velocity (given applied force descriptions). With weight and 
buoyancy acting, we need additional equations (constraints) that will link vehicle attitude 
to motion so that the combined set of equations can be solvable. The constraints to be used 
are developed using the Euler angles which are a general set of angles used to define the 
attitude of the rigid body. From these angles a relationship between rate of change of 
attitude and the body rotational rates given earlier as the angular velocity vector 
(h = [ p,qj J , can be developed. 

8. Incorporation Of Vertical Forces Into The Equation Of Motion 

The weight and buoyant forces that act at the centers of gravity and buoyancy 
must be defined from static analyses. For submerged bodies the weight and buoyancy force 
vectors do not change with vehicle attitude. For a surface ship, B will change with 
attitude, and the righting moments must be computed from Naval Architectural 
considerations. Assuming that weight and buoyancy are fixed in relation to the body fixed 
frame, fF = 0/ + 0 J + (mg)K, and B = 0l + 0J -(pV )K . Since the weight and buoyancy 

terms in the applied forces act in the global vertical direction, they must be transformed 
into components in the vehicle fixed frame before they can be added into the equations of 
motion Returning to equation 2.4, we therefore find the components along the vehicle 
fixed frame to be the third column of that transformation matrix. The net vertical force 
components become. 
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l=(W-B) 



- sin 9 
cos Q sin <}> 
cos 0 cos <J) 



( 2 . 20 ) 



The weight portion of the vertical force acts at the center of gravity of the vehicle. The 
buoyancy portion of the vertical force acts at the center ofbuoyancy. Because these forces 
act at positions away from the body center they exhibit a moment about the body center 
given by, 



m =Wp G x 



-sinB 




- sin 0 


cos Q sin § 


-flp B x 


cos Q sin § 


CO50CO5(f) 




cosQcosfy 



( 2 . 21 ) 



This moment will not be zero even if W and B are identical because it is not usually the 
case that p p and p e are collocated. For static stability it is advisable to locate the center of 
gravity below the center ofbuoyancy The total vertical force vector can be written as. 



Fc- 



fc 



m r . 



( 2 . 22 ) 



and is added negatively to the left hand side of the equations of motion. 

9. Development Of The Full Six Degrees Of Freedom Non-Linear Equations 

Of Motion For A Marine Vehicle 

Define a vector x of vehicle body frame velocities to be, x = [u,v,w,p,q,r ] , 
and a vector z of global positions to be z = [ X ,Y ,Z,<}>,0,V|/ /, then considering M as a 
6x6 mass matrix including translational and rotational inertial elements, the equations of 

motion can be written in the following vector form, 

Mi + f(x)+F g (z) = F h (2.23) 

and, 

:+g(x,z) = 0 (2.24) 
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Therefore with suitable knowledge of the excitation force and moment loads, 
F h , as a function of time and/or vessel motion, a solution for the vehicle's dynamics can be 
obtained. A more detailed insight into the development of the twelve differential equations, 
in first order form given by the foregoing analysis shows, 



m{v + (bxp G } + m((bx(bxp c +6)xvj+ f c ( : ) = 



X, 



(2.25) 



and. 



I 0 {b + m{p G xi’} + (bx(J 0 (b) + m{p G xcox v} + m (:) = 



"f. 



It helps here to define the cross product coefficient matrix so that, 

co x p G = -p G x to = -[ cros(p G )](b 

where. 



(2.26) 



(2.27) 



/cros(p G )7 




(2.28) 



_y G ~ x g 0 . 

Now collecting the mass matrix coefficients into a 6x6 matrix including the inertia cross 
coupling effects, results in the following equation; 
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(2.29) 



The remaining terms on the left hand side of the equations of motion arising from the 
centripetal and coriolis accelerations become, 



/(*) = 



m(ibx(bxp c + (bxv) 
co x (I o (0) + mfp G x w x v} 



(2.30) 



The double vector cross products are nonlinear in the primary velocity variables and hence 
the need for the nonlinear functional, /(•) The reader can perform the indicated 
manipulations to express individual equations within the set if so desired. 

The components of the hydrodynamic and external forces and moments acting 
on the vehicle body are separated in the above analysis into six components each acting 
along the vehicle body fixed coordinate axes and form the total vector of forces and 
moments as, 

F h (t) = [X f (t),Y f (t),Z f (t),K f (t),M f (t),N f (t)]\ (2.31) 



where the vector components in order refer to the surge., sway, heave forces, and the roll, 
pitch, and yaw moments respectively. 

The long form of the six degrees of freedom equations of motion can be written 
as follows, 



m[u-vr + Mq-x G (q z +r) + y G (pq-r) + : c (pr + q)J = 
-(W - BJsin Q+ X f 
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m [ v + ur - wp + x G ( pq+r )- y G ( p 2 + r 2 )+zjqr - p)J = 
(W - B JcosQ sin <{) + Y f 



(2.33) 



m[w-uq + vp + x G (pr-q)+y c (qr + p)-z c (p 2 +q 2 )] = 

(W-B) cosQ costy + Z f ( } 

Kp + (!:- Iy)qr + IJ pr-q)- IJ q 2 -r : )-IJpq + r) + 

m l y g (w-uq + vp ) -: c (v + ur-Mp)J = (y G W - y B B)cosQ cos § (2.35) 

~(: G W -z g B) cos 0 sin (j) + K f 

I y q + ( l x - I.Jpr - IJ qr + p) + IJ pq-r) + IJ p 2 - r 2 ) - 

m[ x G f\i - uq + vp)-z G ( u - vr + wq )]=■-( x G W - x b B)cosQcos§ (2.36) 

~(z c W - z b B ) sin 0 + M f 

I : r + (I y - l z )pq - IJ p 2 -q 2 )- IJ pr + q) + IJ qr-p) + 

m[ x G ( v + ur - wp ) - y G ( u-vr + wp)] - ( x G W- x b BJcosQ sin <{> (2.37) 

+Iy G fV- y B B)sin Q + N f 

while the kinematic relations for the vehicle rate of change of attitude and motion over the 
ocean bottom require equations 2 1 1 and 2.8. 

10. Adaptation Of The Non-Linear Equations Of Motion To The Vertical 
Plane 

Restricting the motions of the vehicle to the vertical (dive) plane, the only 
significant motions that must be incorporated to effectively model the vehicle in the dive 
plane are, the surge velocity («), the heave velocity (w), the pitch velocity ( q ), the pitch 
angle (0) and the global depth position (Z) This restriction simplifies the twelve 
previously developed equations to a system of four non-linear equations of motion, which 
are; 

m( w-uq- z G q- - x G q) = 2 f (2.38) 
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where, 



(2.39) 
(2 40) 
(2.41) 



I y q + m: G wq - mx G (w-uq)= M f 

e = q 

Z --u sin 8 + w co50 



and, 



Z f = Z q q + Z>- + Z q uq + Z„uw - 
1 , ( H ’ -xqf 

,«l I” -^1 



1 ✓ \3 

-p f CJb{x) ^—^ dx + {W-B)cose+u\Z s 8 s + Z B 8 b 

2 \w-xq\ 



(2.42) 



M ^ = M .q + M .w+ M uq + M uw 
f q w q * w 

1 C d 6 ( x ) ^ xdx 

--pJ h’-^l 

L tail 



(2.43) 



-(xJV -x D B)cosQ-{zJV -z B)s\x\ 6+u 2 {M s 8 S + M s 8 b ) 
results from expanding the Z f and the M f terms from (2.31) in a first order Taylor series 

expansion and incorporating both hydrostatic and fluid drag forces. 

Equations 2.38, 2.39, 2.40 and 2.41 can be linearized for a level flight path when 
the dive plane angle is zero, i.e. 5 o =0. By setting all time derivatives to zero, and 
neglecting for the moment the hydrodynamic drag terms the following are obtained; 



Zjiw + (W- B JcosO = 0 (2 44) 

Mjtw - (x Q W - x b B )cosQ-(z g W - : g B)sinQ = 0 (2 45) 

q = 0 (2 46) 

—u sin 0 + h ’ cos 0 = 0. (2 47) 
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If the assumption is made that the vehicle is neutrally buoyant, then it can be said that 
x G = x B and W = B From this equations 2.44 and 2 45 can be reduced to; 



Zjuw = 0 (2.48) 

-: B )BsinQ = 0, (2.49) 

which yield the nominal position w o = q 0 = 0 and sinQ o = 0, resulting in the 0 O solution as 

either 0 O = 0 or 0,, = n. Choosing to linearize around the nominal point 0 O = 0 results in; 

< ? 2 = 2q o q = 0 (2.50) 

wq =M' o q + q o w = 0 (2.51) 

5//»0 = co50 o 0 = 0 (2.52) 

HCO50 = (-M 0 sinQ 0 )Q+(cosf} 0 )w = w . (2.53) 

The linear equations of motion are then written as; 

( m - Z H )\i’ - ( Z q + mx G )q = Zjm + ( Z q + m)uq + Z s u 2 b (2.54) 

( Iy - )q-( M w + mx G )w = M h uw + (M q - mx G )uq -(: G -: B )WQ + M s u 2 b (2.55) 

Q = q (2.56) 

Z = -m0 + w. (2.57) 



As equations 2.42 and 2.43 show, both Z§ and AT§ are a linear combination of 
the respective stem and bow hydrodynamic control surface coefficients and the respective 
input value of 8 This makes the system of equations as a multiple input system. To 
reduce this system into a single input system the linear combination of control inputs will 
be modified into the following form; 

Z b =(Z & + o2 h ). (2.58) 

This will allow a single input 8 to control both stem planes and bow planes, and will 
cause the bow planes to be slaved to the stem planes This technique is known as dual 
control. The value CX will range from -1 to 1. The selection of the value of CX will allow 
the planes to operate as desired for the particular maneuvering condition, i.e., a = 0 for no 
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bow plane control, a = -1 for bow plane and stem plane control opposed to each other, 
yielding the maximum pitch moment, and a = 1 for bow and stem plane control in the 
same direction, yielding the maximum heave force 

In state space form these equations can be written as; 
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where the coefficients a tJ and b t are given by; 

D v =(m- ZJ( 1, - M q ) - ( mx G + Z q )( mx G + M W ) 

a u D v =(I y -M.)Z w + (mx G + Z q )M w 

a ]Z D v = (I y -M q )(m+Z q )+(mx G + Z q )( M q - mx a ) 
a n D v = -( mx G + Z q )W 
b\D v =(I y - M-)Z h + (mx G + Z q )M h 

a 2] D v = (m- ZJM w + (mx G + MJZ w 
a zz D v =(m- ZJ( M q - mx G ) + (mx G + Mj(m + Z q ) 

a zz D v = -(m-Z w )W 
b z D v =(m-ZJM & + (mx G + MJZ h 



and : GB = : G - : B is the metacentric height. 



( 2 . 59 ) 



( 2 . 60 ) 

( 2 . 61 ) 

( 2 . 62 ) 

( 2 . 63 ) 

( 2 . 64 ) 

( 2 . 65 ) 

( 2 . 66 ) 

( 2 . 67 ) 

( 2 . 68 ) 



This linear system of equations becomes the basis for the control law design in 
the following section. 



B. CONTROL LAW DESIGN 
1. Introduction 

The control design problem can be stated as follows: Given the system; 
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x = Ax + Bb, 



(2.69) 



where the state vector equation is. 



0 



U’ 



X = 



9 

Z 



(2.70) 



how do we find 5, such that the system will behave as desired. The type of control that is 
of interest in this problem is closed loop control, where 5 is a function of the state x. 
Since the state x is used to determine the control effort b(x) it is called feedback control. 

2. Pole Placement 

A linear full state feedback control law is introduced in the form 



b = -Kx, (2.71) 

where K is the feedback gain vector to be determined such that the closed loop system of 
equations 2.69 and 2.71 has the desired system dynamics. Substituting equation 2.71 into 
equation 2 69 yields, 

x = ( A- BK )x . (2 72) 

The actual characteristic equation of this closed loop system is given by 

det[A-BK-sl} = 0. (2 73) 

The gain vector K can be chosen such that the actual characteristic equation assumes any 
desired set of eigenvalues. If the desired locations of the closed loop poles are chosen at 
s = s t for /= \, ,n, the desired characteristic equation becomes 

(s-sj(s-s : )...(s-sj = 0 (2.74) 

The required values of K are obtained by matching coefficients in the two polynomials of 
the actual and desired characteristic equations 
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Consider the previously developed linear state matrix equation 2 59. 
Substituting equation 2.71 into equation 2.59 results in the closed loop system 
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(2.75) 



The characteristic equation of the closed loop system is given by 
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a 23 Z GB 



-b^u 2 k x 
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= 0, 



which after some algebra reduces to 

s 4 +(A 2 k 2 + A 2 k 2 - EJs 2 + (-B ] k ] - B 2 k 2 -B 2 k 2 - B A k A -E 2 )sr 
+ (~C ) k ] - C 2 k 2 - C A k A - E 2 )s+'(-D\k A - D 2 k A ) = 0 

where , 

A = -B a =b ] u 2 
Aj = - B , = b 2 u 2 
B 2 = ( a 22 b\ - a u b 2 )u 3 
— Cj — (a^b 2 ci 2 \b\)u 

C 2 D, — (@ 22 b\ — ct\ 2 b 2 )zQgU 
C A — ( a 2 ' P \ "^^2 ~ ^\ 2^2 
A =( a \A ~ a iA) ui 

A =(a u +a 22 )u 

A ^23 “GB "^^^12^21 _ ^11^22 A 



(2.76) 



(2.77) 



(2.78) 

(2.79) 

(2.80) 
(2.81) 
(2.82) 

(2.83) 

(2.84) 
(2 85) 
( 2 . 86 ) 
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^■ 3 — (^n a 2\ < *nCi 2 - i )z CB u (2.87) 

Now if the closed loop poles are placed at -p x , -p 2 , - p 3 , and - p 4 , then the 



desired characteristic equation is, 





(s + p ] )(s + p 2 )(s+p 3 )(s + p 4 ) = 0 . 


(2.88) 


or 




5 4 + a,5 3 + a,5 : + a,5 + a 4 = 0 


(2.89) 


with. 




a \ = P\ + P: + Pi + Pa 


(2.90) 




<*: = PxPz + P\Pi + P\Pa + Pi Pi + P2P4 + P3P4 


(2.91) 




a 3 = P1P2P3 + P1P2P4 + P\PiP* + P1P1P1 


(2.92) 




«4 = PxPlPlPl 


(2.93) 



The control gains can now be computed by equating coefficients of the actual and desired 
characteristic equations, 

AJk 2 + AJc-i = -a, - E x (2.94) 

B x k | + B-M- + B 3 k 3 + B 4 k 4 — ot-, + E 2 (2.95) 

C x k x + C z k 2 + C 4 k A =a i + E i (2.96) 

(D^ + Lhjki =a 4 . (2 97) 

Now that a method for computing the gains of a controllable single input system 
that will place the poles at any desired location has been developed, the selection of pole 
location must be accomplished 
3. Pole Location Selection 

In a typical second order system control law design the transient response 
specifications will be given. This results in an allowable region in the 5-plane where the 
desired location of the poles can be obtained. For higher order systems the concept of 
dominant roots can be employed In selecting poles for a physical system it is necessary to 
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look at the physics of the system If the poles are specified too negative, a very small time 
constant for the control system will result, and the physical system may not be able to 
react that fast. The control law u = -Kx , implies that for a given state x the larger the 
gain, the larger the required control input. In practice there are limits typically placed on u 
(actuator size and saturation). Occasional control saturation is not serious and may even 
be desired. A system that never saturates is in all likelihood over designed. 

Considering the control law design to stabilize the submarine to a level flight 
path at 0 = 0 it will be required that the submarine return to level flight, after a small 
disturbance in 0 or Z, within the time it takes for the vehicle to travel three ship lengths 
Since the submarine is about 14 feet long and it travels at 5 ft/sec, the required recovery 
time is about 10 seconds. This means that the time constant is about 3 seconds, and the 
closed loop poles should be placed at approximately -0.3. 

Control design using pole placement is very easy using MATLAB The 
appropriate MATLAB command is place, which accepts as inputs the A and B matrices 
along with a vector of the desired closed loop poles, and returns the gain vector K. Using 
equation 2.59 with a nominal speed of 5 ft/sec and the vector 

p = [- 0.3 -0.31 -0.32 - 0.337, (place does not like poles in the exact same location), 
the gain vector K was calculated using MATLAB and by simultaneously solving equations 
2.94 through 2.97 yielding the same results. Substituting the gain vector and state variable 
vector into equation 2.71 results in a control law of 

5 = -7-0. 99 1 70 -0.8333m' -0.6026^ + 0.03 5 1Z;. (2.98) 

Using this control law along with equations 2.38 through 2.43 a simulation 
program was developed to investigate the vehicle's response to initial disturbances, the 
results.of which will be presented in the following chapter. 
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ffl. PROBLEM IDENTIFICATION 



A. VEHICLE SIMULATIONS 

With the control law developed and the vertical plane equations of motion defined, a 
vehicle simulation program was developed (see Appendix A) This program was used to 
simulate the vehicle's response at a variety of speeds and initial disturbances. The 
simulations were conducted with several simplifications applied to the vertical plane 
equations of motion (E.O.M.), equations 2.38 through 2.43. The simplifications or 
constraints applied to the E.O.M. were , 

1) to consider the body drag forces negligible, 

2) to consider the vehicle to be neutrally buoyant, i.e., (W=B), 

3) to assume that the locations of the center of gravity and center of buoyancy, 
in the X-Y plane, are coincidental, and 

4) that the rudder is restricted to ± 23 degrees. 

These assumption resulted in the reduction of equations 2.38 through 2.43 to the 
following system of equations; 
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forming the basis of the initial vehicle simulations. 
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The vehicle's response to small disturbances, as stated earlier, was investigated over a 
wide range of speeds. A small sample of these simulations is shown in Figures 3-1 through 
3-3. 

As Figure 3-1 depicts, the vehicle returns to level flight in approximately 10 seconds. 
This response is as expected based on the control system design of Chapter II. The 
simulation results represented in Figure 3-2 also show the vehicle steadying out in a level 
flight at the desired depth, however the vehicle requires a significant amount of time for 
this to occur. The reason for this extended amount of time is that the control law gains 
were set for the nominal operating speed of 5 fps, and the vehicle is being simulated at 
one-third of that speed This results in a very slow return to ordered position caused by 
the reduced hydrodynamic effects at this speed 

The final plot of the sample vehicle simulations, Figure 3-3, illustrates a problem with 
the vehicle stability. As can be seen, the vehicle does not return to the desired position, but 
instead steadies out at a pitch angle of eight degrees, with a steady state dive plane value 
of 23 degrees. This same type of response, but with different pitch angle values, was 
observed in all the simulations conducted below approximately 19 fps These results 
indicate that there is a critical point, or speed at which the vehicle stability is affected 
Determining this point will be the basis for the remainder of this chapter. 

B. CRITICAL SPEED IDENTIFICATION 
1. Eigenvalue Analysis 

Consider the nonlinear system of state equations 3.1, which can be written in the 

form, 

x = f(x). (32) 

It is known that the equilibrium points, x a , of the system are defined by 
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Figure 3-1. Vehicle Response For a Nominal Operating Speed of 5 fps. 



28 



non dimensional depth theta/delta (degrees) 
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Figure 3-2. Vehicle Response at 2.0 fps. 
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Figure 3-3. Vehicle Response at 1.885 fps. 



30 



/(*,) = 0 . 



( 3 . 3 ) 



This is a nonlinear system of algebraic equations and it may have multiple solutions in x 0 , 
which means that the nonlinear system of equations may have more than one position of 
static equilibrium. If one equilibrium value x Q , is selected, the stability properties can be 
established by linearization. The linearization converts equation 3.2 into the form 

x = Ax, (3 4) 

where A is the Jacobian matrix of f(x) evaluated at the equilibrium point x Q , 



-I 



(3.5) 



and the state x has been redefined to designate small deviations from the equilibrium x Q , 

x^x-x a . (3.6) 

As long as all eigenvalues of A have negative real parts, the linear system will be stable. 

This means that the equilibrium x 0 will be stable for the nonlinear system as well. This 
statement is nothing more than Lyapunov's linearization technique. 

The question that must be answered is, what happens if one real eigenvalue of 
the linearized A matrix is zero? The interesting case here is when the rest of the 
eigenvalues have all negative real parts, otherwise x 0 is unstable and the problem is solved. 
If the case of a zero eigenvalue appears too specialized to be of any practical use, consider 
this: Assume that f(x) depends on one physical parameter and that physical parameter is 
allowed to vary over some range. Then clearly A will depend on that parameter and as the 
parameter varies, it is possible that one real eigenvalue of A will become zero for a specific 
value of the parameter. The problem is then to establish the dynamics of the nonlinear 
system as one real eigenvalue of A crosses zero, i.e., goes from negative to positive As 
the solution evolves in time, things are interesting only along the direction of the 
eigenvector that corresponds to the critical eigenvalue. Along the rest of the directions in 
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the state space, everything should converge back to the equilibrium; remember it was 
assumed that all remaining eigenvalues of A have real negative parts. Although, strictly 
speaking, this is a true statement for linear systems, there are technical reasons that force it 
to be true for nonlinear systems as well, the only difference is that the corresponding 
directions in the state space are curved instead of straight. 

By taking the previously developed characteristic equation 2.77, and 
remembering that the roots of this equation are the eigenvalues of interest, it is easy to see 
that if a 4 of equation 2 93, is set equal to zero, then one eigenvalue will be zero. If this is 
done, and recalling that k 4 holds a non-zero value, equation 2.97 reduces to 

^2 1 A “Gb( ^13^2 ^2 1 A ) ~ ^ ' (3 -7) 

Since all parameters of equation 3.7 have a fixed value with the exception of u, there must 
be some value of u that causes the linear A matrix to become unstable. Recognizing this 
fact will allow equation 3.7 to be solved in terms of u yielding 



u = 



: gb (jhP\ a i iPz ) 

( a u b 2 ~a 2] b\) 



-il/2 



(3.8) 



Using equations 2 60 through 2.68, equation 3.8 can be simplified into the following form; 



r -ii/2 

• 7 R 

u= =2^ . (3. 

[(M w Z h -Z w M s \ 

Substituting the appropriate values from Appendix B, equation 3.9 can be solved for the 
value of u that causes the systems to become unstable This results in a value of u = 

1 8979 fps using an a = 0. This value corresponds very closely to the value of 1 .9 fps 
obtained by vehicle simulations 
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As a check to this solution, a MATLAB program was written to compute the 
closed loop eigenvalues of the matrix equation 2.59 as a function of the speed u 
Figure 3-4 shows a plot of the real parts of the eigenvalues versus speed It can be seen 
that the real part of the eigenvalue, represented by the dash-dot line, becomes positive at 
speeds less than approximately 19 fps supporting the earlier findings (For clarity 
purposes, the values of the eigenvalue represented by the dash-dot line were scaled by a 
factor of ten.) 

2. Steady State Analysis 

To determine the existence of multiple steady state solutions, for the system 
represented by equation 2.59, it is necessary to perform a steady state analysis on the 
system of equations which models the vehicle dynamics. Using matrix equation 3.1, and 
setting all time derivatives equal to zero, results in the following set of equations; 



<7 = 0, 


(3 10) 


Z q uq + Zjrw + u 2 Z h h = m(-uq G q 2 ), 


(3.11) 


M q uq + M w uw - Bz gb smd + u 2 M & d = m: G wq , 


(3 12) 


—u sin 0 + w cos 0 = 0 


(3.13) 



Simplifying these four equations yields the following two equations in terms of 0 , w, and 
the physical parameters of the vehicle; 

Z u u 2 tan 0 + w 2 Z 6 5 = 0 (3 14) 

Mji~ tanQ - B: gb sin 0 + u 2 M & h = 0 . (3 15) 

Now if equation 3 . 1 4 is multiplied by M & , equation 3.15 is multiplied by Z 5 and the two 

resulting equations are set equal to each other and simplified, the following single equation 
is obtained; 

[(Z W M S - M h Z & )u 2 + Z & B: gb cosQ ] sin 0 = 0. (3.16) 
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Figure 3-4. Real Part of Closed Loop System Eigenvalues as a Function of Speed. 
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Through inspection it can be seen that equation 3 16 could have multiple 
solutions in 0 , besides the trivial solution of 0 =0 So, if this equation is rearranged and 
solved for 0 *= 0, solutions occur when 



CO.V0 = 






^-8 B-GB 



Therefore, the steady state value of 0 is represented by 



(3.17) 



0, = cos 






B: cb 



(3.18) 



which may have multiple solutions based on the value of the speed u As discussed 
previously, these multiple solutions were evident in the simulations conducted earlier in 
this chapter 

Knowing that the maximum value for cos0 is equal to one, the right-hand side 
of equation 3.17 will be true for values less than or equal to one. If this constraint is 
imposed on 3. 17, then the equation may be manipulated into the following form 



u~ < 






(M W Z S -Z W M S ) 



(3.19) 



Rearranging this equation and solving for the upper limiting case yields the same results as 
equation 3 9 Again this supports the critical speed values discussed earlier 

Equation 3.19 can be expressed in a non-dimensional Froude number form by 
dividing the equation by both the gravitational constant g and the physical parameter : GB 
and taking the square root of the result. This will convert equation 3.19 into the non- 
dimensional form given below; 



Fn = 



jiL< I M 

S~gb ig(M H Z h -Z w MJ 



(3.20) 
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With the steady state 0 solution derived, equations for both the steady state 
value of 5 and steady state value of Z can be determined. Using equation 3.14 and 
substituting the steady state value of 0 will allow the steady state 5 equation to be 
obtained, 

5 „ = (3.2i) 

To achieve the steady state solution of Z, begin by writing equation 2.71 in 
expanded general form as 

5 = -f £,0 + k 2 w + k 3 q + k 4 Z ) (3.22) 

Then by applying the same conditions to equation 3.22 that resulted in equations 3.14 and 
3 15, and using the previously developed steady state equations for 0 and 5, an equation 
of the form 

-7 

— 51 tan 0 M + kfi ss + k z u tan 0^ 

Z a = — (3.23) 

— *4 

is obtained 

The loss of stability of an equilibrium and the generation of additional 
equilibrium states, is called a pitchfork bifurcation and is very common in nature The 
buckling of a beam is one such example Using equations 3.18, 3.21, and 3.23 as a basis 
for another MATLAB program, the steady state solutions of 0, 5, and Z versus Froude 
number were investigated with the results displayed in Figures 3-5 through 3-7 These 
three plots are referred to as supercritical pitchforks, so named because upon the loss of 
stability of the trivial equilibrium the additional nearby equilibrium states are stable 
Graphically, this can be represented in Figures 3-5 through 3-7, where the solid curves 
represent stable and dotted curves represent unstable equilibria [Ref. 17], These three 
figures also demonstrate the control input saturation. Upon control surface saturation, the 
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steady state equations for 8 , 5 and Z are no longer valid due to the limit placed on the 
maximum plane angle The steady state equations that hold for this region of saturation 
are given below; 



= ± 22.9° (0.4 radians ), 



w = 






(3 24) 
(3 25) 



and. 



8„ = sm 



-i 



~GbB 



(3.26) 



As can be seen in Figures 3-5 and 3-6, at an approximate Froude number of 1 05 the 
control surfaces saturate resulting in the linear solutions displayed below that saturation 
Froude number. In Figure 3-7 there is no steady state solution because if the steady state 
values of 5 and 6 are placed into the Z equation, Z is non-zero below the saturation 
Froude number 

A parameter introduced in Chapter II was a, which allowed us to go from a 
multiple input system to a single input system. This parameter has a significant effect on 
the location of the critical speed and/or Froude number Figure 3-8 shows the relationship 
been the critical speed u cr and the metacentric height :q B for different values of a. 

These three curves shown in Figure 3-8 can be converted into a single curve by 
plotting a versus Froude number, and is shown in Figure 3-9. As can be seen in this plot 
the critical Froude number varies from 0.81 to 1.22 over the allowable range of ( X 

In addition to the effect that a has on the location of the critical Froude number, 
it also effects the magnitude of the steady state values at a given Froude number and the 
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Figure 3-5. Steady State Values of 0 vs Froude Number. 
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Figure 3-6. Steady State Values of 5 vs Froude Number. 
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Figure 3-7. Steady State Values of Z vs Froude Number. 
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Figure 3-9. Froude Number as a Function of a. 
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location at which the control surfaces saturate When comparisons are made, it can be 
shown that as the value of a becomes more positive the location of the critical Froude 
number and steady state values of 6 and Z become greater Comparing the changes in the 
5 steady state solutions, the shape of the 5 pitchfork does not change however, the 
location of the critical point moves in the direction of increasing Froude number as a 
increases. The results of this discussion are shown graphically in Figures 3-10 through 
3-12 



An interesting result is displayed in Figure 3-13 This figure depicts the 
relationship between a and the critical or saturation Froude number. At the lower Froude 
numbers, there is very little difference between the critical Froude number and the 
saturation Froude number This demonstrates that upon reaching a low critical Froude 
number the control surface immediately saturates attempting to keep the vehicle stable At 
the higher Froude numbers there is a significant difference between the critical Froude 
number and the saturation Froude number This is because as the Froude number 
approaches 2.78, the value used in the control law design conducted in Chapter II, the 
control surfaces will not have to saturate to maintain stability. 

3. Controllability Analysis 

In general, every system of the form, 

x = Ax + Bu, ~,n\ 

y = Cx, < j27 > 



can be divided through a series of transformations into four subsystems 

1) A controllable and observable part 

2) An uncontrollable and observable part 

3) A controllable and unobservable part. 

4) An uncontrollable and unobservable part 
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Figure 3-10. Comparison of Steady State Values of 6 for Different Values of a. 
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Figure 3-11. Comparison of Steady State Values of 5 for Different Values of a. 



45 



Z (feet) 



STEADY STATE VALUES OF Z vs. FROUDE NUMBER 



50 ' 



40 - 



- for alpha = 0 
-- for alpha = 1 



30 1 



20 f 



10 



for me ta centric height = 0.1 ft 



0L. 



- 10 ' 



- 20 ' 



- 30 ' 



-40 - 



- 50 - 



1.05 1.1 1.15 1.2 1.25 1.3 



Fn 



Figure 3-12. Comparison of Steady State Values of Z for Different Values of a. 
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Figure 3*13. Relationship Between Critical/Saturation Froude Number and (X. 
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This is known as Kalman's decomposition theorem Now, since the transfer function of 
any system is determined by the controllable and observable subsystems only, then the 
transfer function may contain less information than is actually needed to model the 
complete system The precise definition of controllability is: 

A system is said to be controllable if any initial state x(t Q ) can be driven to any 
final state x(tj-) using possibly unbounded control u(t) in finite time t Q < t < tj-. 
From the state equations 3.27, this should depend only on A and B, where A is a n X n 
square matrix The criterion for controllability is as follows: Compute the controllability 
matrix 

c = [B,AB,A-B A n -'BJ, (3.28) 

and the system is controllable if and only if the rank ofc (the number of linearly 
independent rows or columns) is n. Roughly speaking, c shows how possible it is to 
change the state of a system using the input. For a single input system B is n x 1 and c is a 
square matrix The test is then that c be nonsingular 

det(c)* 0. (3 29) 

Now, if the controllability matrix is formed using the A and B matrices from 
equation 2 59 then the following matrix results; 



where, 



C ,1 


C \2 


c o 


C 14 


C 2\ 


C 22 


C 21 


c :4 


C l\ 


C 32 


c 33 


^34 


C M 


C*2 


^43 


C 44 



c„ =0, 

c \z =b 2 u ~ , 

c, 3 = a zl u 3 /), +a 22 u\. 



(3.30) 



(331) 

(3.32) 

(3.33) 
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c, 4 =b i u 2 (a n a^u 2 + a 2l a 22 u 2 ) + b 2 u 2 (a 22 z GB + a 2] a r u 2 + a 22 u 2 ), (3.34) 

c 2] = b x ii 2 , (3.35) 

c ;; = a,,i/Y l + a ]2 u 2 b 2 , (3.36) 

c :3 =6 I i/ : ^a 11 *M : + a ;i a i; Y^ + i ; ;/Ya l3 r cs + a u tf, 2 M : + ar 12 a 22 i / : ), (3.37) 

c ;4 =i | j/Y«n M ^°ii' M: +Y: a :i M: ^ + a 2i w ^ 0 i3 r Gs +a \\ a \z u ‘ + a i2 a 22 u z )) + 

b 2 u z (a n z GB a u u + a ]2 a 2i z GB u + a u u(a u 'u 2 + a i2 a 2] u z ) + (3.38) 

a 22 u(a 22 z GB + a ]2 a 2l u 2 + a 22 u 2 )), 

c 31 =b 2 u z , (3.39) 

c 32 = a 2x u i b x +a 22 ui i b 2 , (3.40) 

c 33 = b } u 2 (a n a 2l u 2 + a 2i a 22 u 2 ) + b 2 u 2 (a 22 z OB +a 2] a i2 u z +a 22 u 2 ), (3 41) 

C 34 b^U ( Cl j 1 2/ ( 47||CJ 2 jW*” "f' Q^CI^U J CJ 2 ^u( ^22 ^ 

b 2 u 2 (a n z GB a 2 ^u + a 22 a 22 z GB u + a n u(a u a 2 \U 2 +a X2 a 22 if )+ (3.42) 

a 22 u(a 2l z GB +a i2 a 2] u 2 +a 22 u 2 )), 

c M - 0, (3 .43) 

c 42 =Z>,m : , (3.44) 

c 43 =a u u 2 b ] +(a r u-u)b 2 u 2 , (3 45) 

and Gy =b ] u 2 (a u : u : +a 2] u(a r u-u)J+b 2 u 2 (a^z GB +a u a r u 2 +a 22 u(a ]: u-u)J . (3.46) 



Now, if the determinant of the controllability matrix is computed as a function of 
Froude number, then it can be shown that the system becomes uncontrollable at the critical 
Froude number The graphical results of this computation are shown in Figure 3-14 for 
different values of the parameter a. 

With the critical parameter that causes the vehicle to lose stability identified, i e , 
the Froude number, and the effect that the parameter a has on the location of this critical 
value understood, the simplifications discussed at the beginning of this chapter can be 
removed from the E.O.M. and an analysis on the effects of incorporating these additional 
parameters can be conducted. This detailed analysis will be the topic of Chapter IV 
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Figure 3-14. System Controllability vs Froude Number. 
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IV. BIFURCATION ANALYSIS 



A. ASYMMETRIC PITCHFORK 

Now that the critical parameter, the Froude number, has been identified, the 
simplifications applied to the E.O.M. in the previous chapter can be removed The 
resulting system of equations is given by equations 2.38 through 2 43 To determine the 
existence of multiple steady state solutions for this new system of equations it is, once 
again, necessary to perform a steady state analysis as conducted in Section B.2 of the 
previous chapter This analysis will reduce the system of equations into the following; 

Z w uw - ^ pC D /lw'|w'| + (W- BJcosB + Z 6 ifb = 0 (4.1) 

Mjiw - ^ pC D x A Aw\w\ - ( x c W - x b B ) cosd ~(: C W - i b B ) sm 0 + M h u 2 5 = 0, (4.2) 

with equation 3.13 also obtained from the steady state Z equation By multiplying equation 
4.1 by M h and equation 4.2 by Z 8 and setting the resulting equations equal to each other 

the following equation is obtained; 



( 4-H - -(~P AC D )(M b - x A Z & )w\w\ + 



(4.3) 



((W - B)M & + (x g W - x b B)Z & ) cos 0 + Z 5 (: G W - z b B) sin 0 = 0 
This equation, along with, 

H’ = ?//a»0, (4.4) 

from equation 3.13 is the exact steady state solution set with all parameters included. 
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Figure 4-1 shows the results obtained when equation 4.3 is solved numerically, using the 
FORTRAN program in Appendix A, for a small value of x GB (the difference in the location 
of the center of gravity and buoyancy), with the submarine neutrally buoyant ( W=B ). 
This small value of x GB perturbs the pitchfork, as seen, from the symmetric pitchfork 
displayed earlier in Figure 3-5. This comparison shows that the critical or bifurcation point 
has moved to a lower Froude number, with the stable solutions represented by the solid 
lines and the unstable solutions represented by the dashed lines in both cases. The stable 
steady state solution, that the vehicle will acquire, is dependent on the initial conditions 
that are given to the vehicle Using the results displayed in Figure 4-1 as an indication that 
x GB may also be a critical parameter, the subsequent task is to determine the relationship 
between the critical Froude number and x GB or f ( Fn,x GB ) . 



B. BIFURCATION GRAPHS 

The first order of business in the attempt to identify f(Fn,x GB ) is to simplify the 
exact pitchfork equation This can be accomplished by replacing the cos0 and s/>?0 terms 
in equation 4.3 with the following Taylor series expansions 



. w 1 u' 3 2 M’tr - 

sin 0 = = 



u 2 u 



2u } 



cose ,,-i^;=5^. 

2 u~ hr 

which will allow the exact solution to be represented as follows; 



(4.5) 
(4 6) 



Z,(: c W-: b B)w z + 2u z C D A(M h - x A ZJw\u] 
-(2u*(ZM h - MJJ + 2irZ & (: G W - : b B))w - hr((W - B)M h , 
+ (x G W - x B B)Z h + ((W- B)M h + (x G W - x b BJZ & )m z = 0 



(4.7) 
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Figure 4-1. Exact Solution Set for Xgb/L=-0.0001. 
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where C D now equals the terms — pC D present in equation 4.3. Defining the following 



parameters 

5 b = B-W, (4.8) 

x G b= x g- x b> (4.9) 

: GB = : G~~B' ( 410 ) 

: a W-:„B = : ot W-z t i, = 5„ (4 11) 

x 0 ir-x B B = x CB H'- x ,5 a =5 I , (4 12) 

X = 1 Fn 2 , (4.13) 

and substituting them into equation 4 7 will allow the coefficients in front of the various 
powers of w to be written as follows; 

H' 3 : A = Z b (z GB W-z B h B ) = Z b (\u 2 W-z B h B ), (4.14) 

u’|h'|: A } =2u 3 C d A(M 6 -x a Z 6 ), (4.15) 

w: A =-2u 2 (u 2 (Z„M 5 -M w Z s ) + Z s (Xu 2 W-z b 5 b )J, (4.16) 

w 0 : A,=-2u 2 (-b B M h -Z b (x GB W-x B b B )) (4.17) 

w 2 :A 0 =(-5 B M b -Z b (x GB fV-x B 5 B J). (4.18) 

The resulting form of equation 4.7 can be written 

A,h’ 3 + /^u’lu’l + Aw + A t + A^m’ 2 =0, (4.19) 



which is a generic pitchfork equation If it is assumed that the submarine is neutrally 
buoyant, and equation 4.19 is divided by Z b z GB W , then the above coefficients can be 

redefined as, 

A 4 = 1, (4.20) 

A =2»C d A(M 5 -x 4 ZJ^-— , (4.21) 

Z J V §“GB 
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(4.22) 






g:, 



GB 



Z & W 



4 = — 2g, 



g= 



GB 



4> = - 



»' X G B g 
g : GB U 



Letting the critical parameter X be defined as 



X = - 



g' 



GB 



will allow equation 4 .19 to be put into the following pitchfork solution form, 

h' 3 + yXh'|m’|-2wY 1 + C^J w ' + 2 Pw 2 A + (3Xw' : = 0, 

where, 



(4 23) 
(4 24) 



(4.25) 



(4 26) 



P = _iss£, (4.27) 

u~ 

r _ g(Z^h ~ M„Z & ) (4 7$) 

s z h w 

y = 2 uC D A(M,-x A Z h )-f- (4.29) 

To determine if the above manipulation and definitions are correct and valid, take the 
symmetric case that was developed in the previous chapter, i.e., x CB =0. This causes 
equation 4.26 to reduce to 

w(w 2 + yX\w\-2u z (\ + t l X)) = 0, (4.30) 

resulting in a solution always occurring at w=0, with two more solutions appearing at the 
bifurcation value 

1 + CX, = 0. (4.31) 

Rearranging and substituting the appropriate values into equation 4.31 will yield 
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ir 



(4.32) 



ZJV 

which is identical to equations 3.9 and 3.20, thus verifying the manipulation 

To find the relationship between Fn and x GB the function h(w,X, ft), and its partial 
derivative hj w,X,P), must be determined. These two functions must then be set equal to 
zero and to each other in order to obtain a function independent of w. As seen, equation 
4 30 is much too complex to complete these requirements, but by neglecting the terms 
Xu vt and Xu ,: in equation 4 26 a simplified pitchfork of the form, 

w i -2u 2 (\ + C,X)w + 2$u 2 X = 0, (4.33) 

or in general terms h(w,X, (5^ = 0, can be obtained which can be manipulated to meet the 
above mentioned requirements, thereby determining the critical (X,$) curve. By taking 

the partial derivative, h M „ of equation 4.33, the following equation results, 

h u = 3w 2 -2u 2 (\ + t,X). (4.34) 

This equation can be set equal to zero, yielding 

w : =|irn + C*J. (4.35) 

which can be substituted into equation 4.33 to obtain the value ofu 1 where the function 
and it's partial derivative is equal to zero; 



3(3X 

w = - . 

2(l + &) 



(4.36) 



This value of u can now be substituted back into equation 4.35 to obtain the desired 
/fX.Py, or 

27P : X : =8i/Vl + CX/. (4.37) 

Attempting to solve this equation analytically for X as a function of (3 , to determine 
the shape of this curve is difficult Therefore, an approximation to this equation is needed 
for small values of P . By rearranging equation 4.37 into the form 
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an equation that determines (3 as a function of X is now produced. If this function is 
represented in a Taylor series expansion and simplified, then equation 4.38 can be written 
as 

*- = \>+ffP 2 ' (“X)) W \ (4-39) 

which can be considered as an approximation to the (X,$) curve, with equation 4.37 as 
the exact (X,$) curve Through inspection it can be seen that the general shape of this 
curve is that of a cusp. 

Figure 4-2 shows a comparison of four bifurcation curves The curves represent the 
relationship between the critical speed U, which can be converted into the non dimensional 
parameter X using equation 4.25, and the value of x GB , which also can be converted into a 
non dimensional parameter (3 using equation 4 27. The curves labeled exact bifurcation 
set, pitchfork bifurcation set, exact cusp and approximate cusp come from the numerical 
solution to equations 4.7, 4.26, 4.37 and 4.39 respectively. Observe that in the general 
vicinity of x GB =0 and near the critical speed the shape of each of these curves is indeed 
that of a cusp, with the peak occurring at the critical speed of 1 9 fps. Now, with the 
comparison conducted in Figure 4-2 complete, equation 4.7, the exact pitchfork solution 
will be used for the remainder of this analysis. 

Returning to equation 4.7, it can be seen that there is only one additional parameter, 
C D , that the cusp curve produced using equation 4.39 does not incorporate (because : GB 
and x GB are incorporated into the non dimensional parameters X and P ). It is therefore 
necessary to conduct a sensitivity analysis to determine the effect that C D has on the shape 
of the cusp curve Figure 4-3 displays the results obtained when this sensitivity analysis 
was conducted. 
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Figure 4-2. Comparison of Bifurcation Curves. 
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Figure 4-3. Comparison of Exact Bifurcation Set for Different Values of Drag 

Coefficient. 
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As seen in the plot, the drag coefficient has no effect on the location of the peak of 
the cusp. The drag coefficient, however, does modify the slope of the cusp as displayed by 
the four cases shown in Figure 4-3. 

C. SOLUTION SETS 

In this section, the equilibrium solutions to equations 4.3 and 4.5 will be investigated. 
As determined in the previous section, the critical parameters that control the location and 
shape of the pitchfork solution sets are X , p and the drag coefficient C D . By running the 
FORTRAN program of Appendix A for several values of x CB L the shape and trend of 
the solution sets can be determined. The result of this comparison, for a C D = 0.0, is 
displayed in Figure 4-4. This plot supports the symmetry about x CB = 0 displayed in the 
earlier cusp curves Observe that the critical or bifurcation point moves to a lower value of 
X as x GB increases, and that for low values of X the equilibrium solutions converge to the 
same value independent of the value of x GB . 

When the comparison is made regarding the effect that the drag coefficient has on the 
solution sets, for a given value of x GB , Figure 4-5 is produced. This variation in drag 
coefficient for a constant value of x GB has a similar effect on the movement of the 
bifurcation point as fixing the drag coefficient and varying x GB . The difference is, however, 
evident by studying the stable solutions The stable solutions for non-zero drag 
coefficients become more linear as the value of C D increases. This results from the 
increased effect of the w|w| term in equation 4 3 as the value of C D increases. 

Although Figure 4-5 shows the actual solution to the exact pitchfork equation it is 
not realistic since the effect of dive plane saturation is not considered. Once the dive 
planes saturate, the equations used to obtain the steady state solutions displayed in 
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Figure 4-4. Comparison of Exact Solution Set for Different Values of Xgb/L. 
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Figure 4-5. Comparison of Exact Solution Set for Different Values of C D . 
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Figures 4-4 and 4-5 are no longer valid Instead, equations 4.1 and 4 2 are modified, still 
remembering that the assumption of neutral buoyancy remains in affect, resulting in the 
following two equations, 

Zjm’ - C d Aw |w| + u~Z b 5 MI = 0, (4.40) 

Mjuw - C d x a Aw\w\ - Wx GB cosQ- : GB W sinQ + M & irb sal = 0. (4.41) 

To obtain an equation that can be used to compute the steady state 0 solution while 
taking into consideration dive plane saturation requires solving equation 4.40 for w and 
then substituting this value into equation 4.41 which should yield an equation independent 
of vi’ The ability to perform this manipulation is hindered by the inclusion of the absolute 
value term in both of the above equations. However, if the absolute value terms are 
neglected by considering the case of C D = 0 . 0 , then a single equation for the steady state 
value of 0 may be obtained. 

Neglecting the drag affects, and performing the algebra will yield an equation of the 
following form, 

( K z & ~ M h ZJir 6 m , + ZW(x GB cosQ + z GB 5//; 0) = O, (4.42) 

which may be used to determine the values of 0 once the dive planes have saturated The 
validity of neglecting the drag terms is questionable, and an analysis that incorporates C D 
will be conducted later in this section once the effect of varying x GB is determined Figure 
4-6 displays the changes that occur to the solution set of Figure 4-1 when saturation is 
considered. Notice that saturation occurs when the solution set is at approximately nine 
degrees, and that when multiple solutions occur, the stable equilibrium states are 
represented by the saturation portions of the solution set. 
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Figure 4-6. Exact Solution Set for Xgb/L=-0.0001 With Saturation Included. 
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When the saturation condition is applied to the previously displayed solution sets for 
different values of x GB (Figure 4-4), the results shown in Figure 4-7 are produced Items 
of interest in Figure 4-7 are that: 

1) The bifurcation point on each solution set moves to a lower value of X as 
the absolute value of x GB increases. 

2) The range where multiple solutions exist decreases as the absolute value of 
x GB increases 

3) The amount of trim produced, in still water, by the given value of x GB can be 
obtained from the point where the upper an lower solution branches 
converge. 

Incorporating the saturation condition into the solution sets raises some additional 
concerns. To begin, what if the point of dive plane saturation occurs after the bifurcation 
points identified in Figure 4-4? If this condition did occur, referring to Figure 4-1 which 
displays the stable and unstable solution branches, there would be at least one value of X 
where there were three stable solution branches. To determine if this condition is possible 
a comparison of solution sets with and without saturation is necessary. This comparison is 
shown in Figure 4-8. As this plot depicts, the range between the bifurcation point and 
saturation point increases as the value of x GB increases. This result will allow the fore 
mentioned concern, of three stable solution branches, to be disregarded. 

A second concern is the validity of the previously developed bifurcation cusp curves 
Since the exact bifurcation equation, equation 4-7, does not include saturation, the cusp 
curves produced using this equation will not accurately predict the occurrence of multiple 
solutions However, if equation 4.42 is expanded in a Taylor series, it can be manipulated 
into a form that will allow a cusp curve, which includes saturation, to be developed. 
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Figure 4-8. Comparison of Bifurcation Points and Dive Plane Saturation Points. 
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Performing this manipulation will produce the following equation; 



if = 



-XofJV + ZoeWZ ,S,„. 



(Z.M.-M.ZJ(&/\ + (Z, zj&r-y 



(4.43) 



By converting if and x GB into the non dimensional parameters X and P the cusp 
curve shown in Figure 4-9 can be produced This cusp displays the saturation point for 
both solution branches, and as seen it differs quite extensively from the exact bifurcation 
cusp without saturation If Figure 4-9 is replotted to allow the peak area of the cusp to be 
more accurately represented, then Figure 4-10 is produced. By using these two cusp 
curv es then the general form of the resulting solution set for any path through the cusp can 
be predicted It is this prediction that the subsequent section will address. 



D. PATH FORMULATION 

The ability to accurately predict how the vessel will respond to various changes in 
operating conditions is critical for proper control of the vessel. It is this prediction that 
necessitates the need for path analysis to be conducted If the ballast condition on the 
vessel is held constant, x GB =-0.0001, and the speed of the vessel is varied, path A of 
Figure 4-11 is obtained. Neglecting for the moment the saturation cusp shown in 
Figure 4-11, it would be expected that a single steady state solution would exist for all 
values, of the parameter sqrt(X), greater than the point where the path intersects the solid 
cusp curve At the speed where the path intersects the cusp the bifurcation point occurs, 
and for speed values which place the path inside the cusp multiple steady state solutions 
occur The results of this path analysis are depicted in Figure 4-1, presented earlier in this 
chapter 

If attention is returned to the path through the saturation cusp of Figure 4-11, then 
the same type of predictions made for the path through the solid cusp can be made for the 
this cusp, with some slight differences. As the speed decreases the path intersects the 
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Figure 4-11. Path Formulation for Xgb/L=-0.0001. 
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upper portion of the dotted cusp at an approximate value of 1.08. At this value saturation 
occurs on one of the solution branches but there remains only one steady state solution 
since the path has not yet entered the cusp As the speed continues to decrease, the path 
intersects and enters the saturation cusp at an approximate value of 0 99 with saturation of 
the other solution branch and the existence of multiple steady state solutions occurring. 
These results are evident in Figure 4-6 

Now, if the speed of the vessel is held constant, and the loading conditions on the hull 
are changed, then the paths shown in Figure 4-12 are obtained Path B, at a Froude 
Number of 1.1, is outside the cusp, therefore only a single steady state solution for each 
value of x GB L should exist While path C, at a Froude Number of 1.0, passes through the 
cusp at ± 0.0001 values oix GB L resulting in single steady state solutions for values less 
than or greater than ± 0.0001, while values between ± 0.0001 result in multiple steady 
state solutions. Through numerical computation the solution set as a function of x GB L is 
obtained with the results shown in Figure 4-13 

As seen in Figure 4-13, path B, represented by the dotted line, does in fact result in a 
single steady state solution throughout the range of x GB -L, however, for path C this is not 
the case. Path C, represented by the solid line, exhibits multiple solutions between the 
x GB L values of ± 0.0001, as predicted. This is a hysteresis curve, with the unstable 
solutions occurring between the comers of the curve. If the value of x GB L began 
at -0.0005 and progressed to the positive value of 0.0001, the solution set would remain 
on the curve with 0 taking on the values shown Once the value of x GB L becomes greater 
than 0 0001 the solution jumps from -9° to 7.5° This same jump would be evident if the 
value of x GB L began positive and progressed negative. These jumps indicate that if the 
dive planes are held constant, then the vessel will instantaneously change its pitch angle. 
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Figure 4-12. Path Formulation at a Constant Speed With Xgb/L Varying. 
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To counter this effect the dive plane angle must be reversed It is this dive plane reversal 
that the following chapter will discuss 

Before proceeding to the next chapter it is necessary to determine the effect that drag 
has on the previously developed bifurcation cusps Using equations 4 40 and 4.41, and 

substituting equation 4.4 into both equations will result in the following two equations; 

Zji 2 tan0-C /y 4w : tan0jtan6) + Z i5 Jr :! 5 Mr = 0, (4 44) 

Mjiw-C d x a Au 2 tan 6jtan 0j - Wx GB cosd- z CB W sin 0+ M s u 2 8 sal = 0. (4.45) 

By using the MATLAB function J. zeros contained in the program of Appendix A, the 
value of d that solves equation 4 44 can be determined as a function of speed. This 
relationship can then be used to solve equation 4 45 to obtain the relationship between x GB 
and u 

With this relationship obtained the cusp curves of Figure 4-14 can be plotted It can 
be seen that for a given value of x GB L increasing the drag coefficient forces the 
bifurcation point to a lower speed. An interesting item to note is the diamond shape peak 
on these cusps This is a result of saturation not occurring, for a small range of speeds, in 
the small region around x GB /L=0. Path formulation through these cusps can be developed 
as shown earlier in this section. 
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Figure 4-14. Bifurcation Cusp With Saturation Included for Different Values of C D . 
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V. DIVE PLANE REVERSAL 



A. STERN PLANE REVERSAL 

The ability of a submarine to avoid detection and perform its mission is dependent on 
its ability to maintain ordered depth. Modem asymmetric submarines exhibit an inherent 
phenomenon known as dive plane reversal, which is difficult to deal with. To introduce 
this physical occurrence, refer to Figure 5-1 [Ref. 14] and consider the case of wanting to 
dive the vehicle to a deeper depth. At moderate to high speeds the stem planes would be 
deflected by an angle 5 S , this would in turn cause a force and moment to be developed due 
to the angle of attack of the planes. The vessel will respond by pitching downward, and 
will assume some angle of attack to the fluid flow causing hull forces and moments to 
develop. The stem plane and hull forces bring the ship to some steady state pitch angle 
which is obtained when the restoring moment and pitching moments balance The vessel 
now has the ability to fly itself to the ordered depth with the normal hull force assisting in 
pushing the vessel to a deeper depth. 

At low speed we would expect the same to be true, however it is not. If the stem 
planes are set to an angle 8 S , this would cause a force and moment to be developed due to 
the angle of attack of the planes, the vessel will pitch downward creating a pressure field 
around the vessel due to the angle of attack with the fluid resulting in hull forces and 
moments to develop. The vessel once again steadies out at some steady state pitch angle 
when the restoring moment and pitching moments balance. Since the speed of the vessel is 
not as great as in the previous case, the hull forces and moments are not as large, therefore 



77 



r * normal force due to angle oe attack on hull 

Zr t NORMAL FORCE DUE to STERNPLANE DEFLECTION ANGLE (6il 

b , * 

M W PITCH MOMENT DUE TO ANGLE OF ATTACK ON HULL 

Mr 6, PITCH MOMENT DUE TO STERNPLANE DEFLECTION ANGLE (&i) 

0 \ 

M^ METACENTRIC PITCH MOMENT DUE TO PITCH ANGLE (0) 

V 





MODERATE and HIGH SHIP SPEED 



Figure 5-1. Stern plane Reversal . 



instead of the vessel flying itself to the ordered depth with the hull force assisting, the 
stem plane force actual pushes the vessel upward causing it to rise instead of dive, with 
some non-zero pitch angle. [Ref 14] This condition was evident in the simulation shown 
in Figure 3.3. 

With this description of the vehicle response complete, it is obvious that to dive the 
vessel at low speeds, the opposite stem plane angle must be ordered This will result in the 
vessel being pushed to the ordered depth at some bow up attitude This required shift in 
dive plane angle is what is known as STERN PLANE REVERSAL. 

To tie together the previously completed bifurcation analysis with the physical 
occurrence just discussed. Figure 5-2 is produced. Figure 5-2 presents the pitch angle and 
stem plane angle required for straight and level flight with the bow planes centered As the 
vessel's speed decreases, the metacentric moment ( Mjd ) makes it increasingly difficult to 

maintain a positive pitch angle on the hull. The stem plane angle is negative trying to 
produce a pitch-up moment, but the downward force of the stem planes requires an even 
larger pitch angle [Ref. 14], Between the speeds of 1.08 and 1.18 knots, the required stem 
plane angle is beyond the normal operating range and the vessel cannot be flown straight 
and level. Below the critical speed, the stem plane and vessel pitch angle must reverse sign 
to maintain neutral trim. 

In a final effort to show the physical significance of stem plane reversal, the force 
input to the stem planes is plotted If the control gains are allowed to change as the speed 
of the vessel changes. Figure 5-3 is produced. This plot shows the non-dimensional force 
per unit depth error input to the stem plane as a function of speed It is basically the 
control system gain used to set the stem plane angle. As can be seen, as the critical speed 
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Figure 5-2. Variation of Pitch and Stern plane Angles as a Function of Speed. 
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Figure 5-3. Control System Force Input to the Stern planes. 
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is approached the value of input tends to infinity requiring the stem planes to operate 
outside normal values. As the speed passes the critical value the sign of this input changes 
requiring the stem planes to operate in the opposite direction. This reversal in sign is 
evident of the stem plane reversal phenomenon shown and discussed earlier in this 
chapter. 

B. BOW PLANE REVERSAL 

The reversal effect that was shown to occur with the stem planes can also be shown 
to occur with the bow planes. If the value of a is allowed to tend to oo, indicating no stem 
plane control and only bow plane control, then the same analysis conducted throughout 
this thesis for stem plane reversal can be performed for bow plane reversal. If the critical 
Froude number, (sqrt(X)), is plotted as a function of a, then the results shown in 
Figure 5-4 are produced This plot shows that for bow plane control only, the critical 
value approaches 2.05 as a approaches oo. 

To show that the same type of results can be obtained for the bow plane analysis, the 
case of neutral buoyancy, x GB = 0 and dive plane saturation, with bow plane control only 
was investigated Referring to Figure 5-5, it can be seen that at a Froude number of 2.05 
multiple steady state solutions occur. These results are similar to the results displayed in 
Figure 3-5 The major differences are that the critical point occurs at a Froude number of 
2 05 vice 1.05, and the maximum value of 0 is only ± 4.5° vice ± 9.5°. The reason for 
these changes is due to the values of Z Sh and M St being only one-half the respective value 

for the stern planes 
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Figure 5-4. Identification of Critical Froude Number for Bow plane Control. 
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C. BUS EFFECTS 

As a final area of analysis, the condition of operating near the surface, such as during 
periscope operations, must be investigated to determine the effects this will have on the 
previously developed bifurcation cusps Operating near the surface will require the 
incorporation of a suction force and moment into the equations of motion. This force and 
moment is a function of the distance away from the surface, and the angle the vessel 
makes with the surface. Therefore, as a simple model of suction effects consider the case 
of operating with excess buoyancy. To conduct the analysis requires the use of equations 
4.1, 4.2 and 4.4. Performing the same substitution and solution technique as was used in 
Chapter V to obtain the saturation curves that incorporated the drag coefficient will result 
in solving the following equations to obtain the bifurcation cusps that would quasi-model 
surface effects; 

Zjiw - - j pCpA u 2 tan 0j tan 0| + ( W - B ) cos 6 + Z s u 2 5 sal = 0 , (5.1) 

A ijiw — -pC D x A Au 2 tan Sjtan 0j - (x G W - x b B) cos 6 

2 (5.2) 

~(z g W-z b B) sind+Mju 2 8„=0 

If the value of C D is set equal to zero then the saturation cusp produced in 
Figure 4-10 can be compared to the results obtained for several cases of excess buoyancy 
Since it is the saturation cusp that controls the location of the bifurcation point, as shown 
in previous the chapter, the non saturation cusp will be neglected for the remainder of the 
analysis. Figure 5-6 shows the results obtained as excess buoyancy is increased As can be 
seen, the cusp peak moves to down and to the right occurring at a lower speed and at a 
non zero value of x GB L. The body of the cusp also shifts to the right, as seen, with the 
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Figure 5-6. Bias Effects Caused by Operating Near the Surface. 
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entire cusp occurring in the positive region of the x GB L range for a value of 0.05 % 
excess buoyancy. 

To determine how the value of drag coefficient effects the cusp curves, consider the 
case of C D = 0.3 shown in Figure 5-7. These cusps were produced for the same values of 
excess buoyancy used in Figure 5-6. Notice that the cusp curves shift to the right by the 
same value seen in Figure 5-6, but the peak value now occurs at an even a lower speed 

In order to observe the sensitivity of the biased cusp to drag only, consider the case 
of setting the excess buoyancy to 0.01 % and varying the drag coefficient shown 
Figure 5-8 These results indicate that increased drag causes the critical speed to occur at 
a lower value They also show that drag tends to counter the biasing effect that excess 
buoyancy was previously show to have. This is evident by the peak point drifting to the 
left as the value of drag coefficient increases. 

In summary, it can be stated that operating near the surface will tend to reduce the 
critical speed determined in previous chapters and the vessel will have to be placed in a 
different trim condition to counter surface effects, and that increasing drag also decreases 
the critical speed but tends to counter the trim condition required. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

• A method for analyzing submarine depth keeping and dive plane reversal has been 
developed and presented. This analysis identified the three phenomena, a positive real 
eigenvalue, a steady state pitchfork bifurcation, and an uncontrollable system, that lead to 
vessel instability 

• The critical parameters have been identified that control the bifurcation. These 
parameters combine speed, loading conditions and hydrodynamic coefficients into 
non-dimensional values that may be applied to other models or vessels 

• Since the non-dimensional Froude numbers are based on the metacentric height 
and LCG-LCB separation, Froude scaling can be used to apply this analysis to full scale 
vessels, allowing for the reduction of model testing during design phases. 

B. RECOMMENDATIONS 

Some recommendations for further research are as follows 

• Develop and incorporate into the bifurcation analysis a more complete expression 
for the free-surface effects. 

• Conduct a Hopf bifurcation analysis to identify the critical parameters that lead to 
vehicle instability at high speeds 
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APPENDIX A 



COMPUTER PROGRAMS 

% THIS PROGRAM IS THESIS l.M. IT SIMULATES THE RESPONSE OF THE 
DTRC SUBOFF 

% MODEL. IT PLACES THE POLES, CALCULATES THE CRITICAL SPEED, AND 
USES 

% PROGRAM THESISDE.M TO SIMULATE THE VERTICAL PLANE REPONSE. 



clg; 

global zg m zgb u U Zq Zw Zdlt Mq Mw B 1 Mdlt K 1 M bm xL Zval Mval rho 

% ESTABLISH GEOMETRIC PARAMETERS 
W= 1 556. 2363 ,B 1 = 1 556. 2363, 

Iy=56 1.32; 
g=32.2; 
m=W/g, 
rho=1.94, 

L=13.9792; 
xg=0;zg=0. 1; 
zgb=0. 1, 

% BEAM AND LENGTH DEFINITION FOR TRAPAZIODAL RULE USE 
bm=[0 .485 .658 .778 .871 .945 1.01 1.06 1.18 1.41 1.57 1.66 1.67 1.67 1.67 1.63... 
1.37 .919 .448 .195 .188 .168 .132 .053 0]; 

xl=[0 .1 .2 .3 .4 .5 .6 .7 1 2 3 4 7.7143 10 15.1429 16 17 18 19 20 20.1... 

20.2 20.3 20.4 20.4167]; 
xl=(L/20)*xl; 
xL=xl-L/2; 

% FACTORS TO CONVERT FROM NON-DIMENSIONAL VALUES TO 

DIMENSIONAL VALUES 

ndl=.5*rho*L A 2; 

nd2=.5*rho*L A 3; 

nd3=.5*rho*L A 4, 

nd4=.5*rho*L A 5; 



91 



Alpha=0; 



% NON-DIMENSIONAL HYDRODYNAMIC COEFFICIENTS 
Zqdnd=-6.33e-4,Zwdnd=-l 4529e-2,Zqnd=7. 545e-3;Zwnd=- 1.391 e-2; 
Zds=(-5.603e-3);Zdb=0.5*(-5.603e-3);Zdltnd=(Zds+Alpha*Zdb); 

Mqdnd=-8 8e-4;Mwdnd=-5.61e-4;Mqnd=-3 702e-3;Mwnd=l 0324e-2; 
Mds=(-0.002409),Mdb=0.5*(0.002409);Mdltnd=(Mds+Alpha*Mdb); 

% CONVERTION OF COEFFICENTS INTO DIMENSIONAL VALUES 

Zqd=nd3*Zqdnd,Zwd=nd2*Zwdnd;Zq=nd2*Zqnd;Zw=ndl*Zwnd, 

Zdlt=ndl*Zdltnd, 

Mqd=nd4*Mqdnd;Mwd=nd3*Mwdnd;Mq=nd3*Mqnd;Mw=nd2*Mwnd; 

Mdlt=nd2*Mdltnd; 

% FORMULATION OF A AND B MATRIX ELEMENTS 

Dv=(m-Zwd)*(Iy-Mqd)-(m*xg+Zqd)*(m*xg+Mwd); 

al lDv=(Iy-Mqd)*Zw+(m*xg+Zqd)*Mw; 

al2Dv=(Iy-Mqd)*(m+Zq)+(m*xg+Zqd)*(Mq-m*xg), 

al3Dv=-(m*xg+Zqd)*W; 

b 1 Dv=(Iy-Mqd)*Zdlt+(m*xg+Zqd)*Mdlt; 

a2 1 Dv'=(m-Zwd)*Mw+(m*xg+Mwd)*Zw; 

a22Dv=(m-Zwd)*(Mq-m*xg)+(m*xg+Mwd)’ , ‘(m+Zq); 

a23Dv=-(m-Zwd)*W; 

b2Dv=(m-Zwd) , ‘Mdlt+(m*xg+Mwd)*Zdlt; 

all =a 1 lDv/Dv;al2=al2Dv/Dv;al3=al3Dv/Dv; 
a2 1 =a2 1 Dv/Dv;a22=a22Dv/Dv;a23=a23Dv/Dv; 
b 1 =b 1 Dv/Dv;b2=b2Dv/Dv; 

% ESTABLISHMENT OF POLE LOCATIONS 
p 1 =[-. 3 - 31 -.32 -.33]; 

% CALCULATION OF CRITICAL SPEED 
Ucr=sqrt(-(al3*b2-a23*bl)*zgb/(al 1 *b2-a21 *bl)) 

% INPUT OF SPEED WHICH THE VESSEL WILL USE IN THE SIMULATION 
U=input('enter speed of vehicle ') 

% CALCULATION OF FROUDE NUMBER 
Fn=sqrt(U A 2/(zgb*g)) 



92 



% CALCULATION OF MATRICES FOR POLE PLACEMENT AND 

CONTROLABILLITY ANALYSIS 

u=5; 

a=[0 0 1 0;al3*zgb al l*u al2*u 0;a23*zgb a21*u a22*u 0;... 

-u 1 0 0]; 

al=[0 0 1 0;al3*zgb all*Ucr al2*Ucr 0;a23*zgb a21 *Ucr a22*Ucr 0, 
-Ucr 1 0 0]; 



b=[0;bl *u A 2;b2*u A 2;0]; 
bl=[0;bl*Ucr A 2;b2*Ucr A 2;0], 

Kl=place(a,b,pl); 

C=rank(ctrb(a 1 ,b 1 )), 

H=(Mw*Zdlt-Zw*Mdlt)/(zgb*B 1 *Zdlt); 
thetal=acos(H*U A 2)* 180/pi; 

Dss=-(Zw/Zdlt)*tan(acos(H*U A 2))* 1 80/pi; 
Z=(Dss+Kl(l)*acos(H*U A 2)+Kl(2)*U*tan(acos(H*U A 2)))/ 

(-Kl(4)); 

% ESTABISHES VEHICLE MASS MATRIX 

M=[l 0 0 0;0 (m-Zwd) -Zqd 0 ;0 -Mwd (Iy-Mqd) 0 ;0 0 0 1]; 

% SET TIME LIMITS AND INITIAL CONDITIONS FOR ODE SOLUTION 
t0=0;tfl =2500*L/U; 
x0=[0 0 0 1]; 

% SOLVES VEHICLE SIMULATION ODE'S 
[t 1 ,X 1 ]=ode45('thesisde',t0,tfl ,x0); 

% REDEFINE THE OUTPUT VECTOR ELEMENTS AND FORM DIVE PLANE 
RESPONSE 

wl=Xl(:,2);ql=Xl(:,3);zl=Xl(:,4);thl=Xl(:,l), 

d=-(Kl(l)*thl+Kl(2)*wl+Kl(3)*ql+Kl(4)*zl); 

ll=length(d); 

for n= 1:11, 
if d(n)> 4, 
d(n)=4, 
elseif d(n)<-.4, 
d(n)=-.4; 
end, 
end; 
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% PLOT THE RESULTS 
subplot(2 11), 

plot(t 1 *U/L,th 1 * 1 80/pi,t 1 *U/L,d* 1 80/pi,'--');title('THET A/DELTA vs. TIME '); 
xlabel('non-dimensional time ');ylabel('theta/delta (degrees)'), grid, 
gtext('— delta');gtext(['for a speed of' num2str(U) ' fps']); 

plotftl *U/L.zl/L);title('DEPTH vs TIME'),xlabel('non dimensional time'); 

ylabei('non dimensional depth ');grid, 

pause, 

subplot(lll), 

% THIS IS PROGRAM THESISDE.M. IT CONTAINS THE SUBOFF MODEL 
SYSTEM OF 

% EQUATIONS TO BE SOLVED BY PROGRAM THESIS l.M 

% SET SYSTEM OF EQUATIONS FUNCTION 
function Xldot=thesisde(tl,Xl), 

% SET CONTROL LAW 

dl=-Kl(l)*Xl(l)-Kl(2)*Xl(2)-Kl(3)*Xl(3)-Kl(4)*Xl(4); 

% APPLY DIVE PLANE SATURATION 
if dl> 4, 
dl=4, 

elseif dl<- 4, 
dl=- 4, 
end, 

% INCLUDE DRAG TERMS AND CALCULATE AREA AND CENTROID 
CD=0; 



if X 1 (2)=0, 

Xl(2)=le-5, 

end; 

if Xl(3)=0, 

Xl(3)=le-5; 

end, 

Zval=bm *((X 1 (2)-xL. *X 1 (3)) A 3)./(abs(X 1 (2)-xL. *X 1 (3))); 
Mval=bm.*(((Xl(2)-xL.*Xl(3)). A 3).*xL)./(abs(Xl(2)-xL *X1(3))); 



ZDragval=0;MDragval=0; 
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% trapaziodal integration 
for n= 1 :length(xL)- 1 , 



Zdragval=0. 5 * (Zval(n)+Zval(n+ 1 ))*(xL( n+ 1 )-xL(n)); 

Mdragval=0.5*(Mval(n)+MvaI(n+l))*(xL(n+l)-xL(n)); 

ZDragval=ZDragval+Zdragval, 

MDragval =MDragval +M d ragval ; 

end; 



Zdr=(0.5)*rho*CD*ZDragval, 

Mdr=(0.5)*rho*CD*MDragval; 

% VEHICLE SYSTEM OF EQUATIONS 

F1=[X1(3); m*(U*Xl(3)+zg*Xl(3) A 2)+Zq*U*Xl(3)+Zw*U*Xl(2)+U A 2*Zdlt*dl-Zdr; 

m*(-zg*Xl(2)*Xl(3))+Mq*U*Xl(3)+Mw*U*Xl(2)- 
B1 *zgb*sin(Xl(l))+U A 2*Mdlt*dl-Mdr; 

-U*sin(X 1 ( 1 ))+X 1 (2)*cos(X 1(1))]; 



mpi=inv(M); 

Xldot=mpr*Fl; 
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c PROGRAM STEADY 
C 

C STEADY STATE SOLUTIONS OF SUBMARINES 
C IN THE DIVE PLANE AT LOW SPEEDS 
C 

C INPUT DATA 
C 

C ICON = 1 : COMPUTATION OF SOLUTION SETS (S.S) 

C 2 : COMPUTATION OF BIFURCATION GRAPHS (B.G ) 

C 

C IVAR = 1 : Fn VARIATION 
C 2 : Xgb VARIATION 

C 

REAL MW, LENGTH,MASS,MDS,MDB,MD, LAMBDA 
C 

DIMENSION B(25),X(25),BX(25) 

C 

C GEOMETRIC PROPERTIES AND HYDRODYNAMIC COEFFICIENTS 
C 

RHO = 1 94 
LENGTH= 13 9792 
WRITE (*,*) ' ENTER CDZ* 

READ (V)CDZ 
CDZ = CDZ*0.5*RHO 
BUO= 1556 2363 

WRITE(V) 'ENTER VALUE OF EXCESS BUO' 

READ (*,*) FAC 
WEIGHT=FAC*BUO 
MASS = WEIGHT/32.2 
XB =0.0 
ZB =0.0 

WRITE (*,*) ’ENTER THE MET ACENTRIC HEIGHT ZGB' 

READ (*,*) ZG 

ZW =-0.013910*0 5*RHO*LENGTH* *2 
ZDS =-0.005603*0. 5*RHO*LENGTH* *2 
ZDB =-0.005603* 0.25 *RHO*LENGTH* *2 
MW = 0 010324*0. 5*RHO*LENGTH**3 
MDS =-0.002409*0. 5 *RHO*LENGTH* *3 
MDB = 0.002409*0 25*RHO*LENGTH**3 
C 

WRITE (*,*) 'ENTER THE VALUE OF ALPHA' 

READ (*,*) ALPHA 
ZD=ZDS+ALPHA*ZDB 



96 



MD=MDS+ALPHA*MDB 

C 

C DEFINE THE LENGTH X AND BREADTH B TERMS FOR THE 

INTEGRATION 

C 

X( 1)=20.4167 
X( 2)=20.4000 
X( 3)=20.3 
X( 4)=20.2 
X( 5)=20. 1 
X( 6)=20.0 
X( 7)= 19.0 
X( 8)= 18.0 
X( 9)=17.0 
X(10)=16 0 
X(1 1)=15 1429 
X(12)=10.0 
X( 1 3 )=7. 7 1 43 
X(14)=4.0 
X(15)=3.0 
X(16)=2.0 
X(17)=l 0 
X(18)=0.7 
X(19)=0 6 • 

X(20)=0.5 
X(2 1)=0.4 
X(22)=0.3 
X(23)=0.2 
X(24)=0 1 
X(25)=0.0 
C 

B( 0=0.00000 
B( 2)=0.03 178 
B( 3)=0. 07920 
B( 4)=0. 10074 
B( 5)=0. 1 1243 
B( 6)=0.1 1724 
B( 7)=0. 26835 
B( 8)=0. 55025 
B( 9)=0.8 1910 
B(10)=0. 97598 
B(1 1)=1 00000 
B( 12)= 1.00000 
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B(13)=l. 00000 
B(14)=0. 99282 
B(15)=0 94066 
B(16)=0 84713 
B(17)=0 70744 
B(18)=0.635 14 
B(19)=0 60352 
B(20)=0. 56627 
B(21)=0. 52147 
B(22)=0 46600 
B(23)=0. 39396 
B(24)=0. 29058 
B(25)=0. 00000 
C 

DO 1 1=1,25 

B(I)=B(I)* 1.6667 

X(I)=(-X(I)+10.0)*LENGTH/20.0 

BX(I)=B(I)*X(I) 

1 CONTINUE 

CALL TRAP( 1 8,B,X,AREA) 

CALL TRAP( 1 8,BX,X,CH) 

XA=CH/AREA 

C WRITE(*,*) 'AREA EQUALS ',AREA 
C WRITE(V) 'XA EQUALS ',XA 
IFC=0 
C 

C OPEN RESULTS FILES 

OPEN ( 1 0,FILE='R 1 O RES', ST ATUS= r NEW') 
OPEN (1 1,FILE-R1 l.RES',STATUS='NEW') 
OPEN ( 1 2,FILE='R 1 2.RES',STATUS='NEW') 
OPEN ( 1 3,FILE='R1 3.RES',STATUS= r NEW') 
OPEN (14,FILE='R14.RES',STATUS =r NEW') 
OPEN ( 1 5,FILE='R1 5.RES , ,STATUS= r NEW l ) 
OPEN ( 1 6, FILE- R 1 6.RES',STATUS= r NEW l ) 
OPEN (20, FILE- C20. RES', STATUS='NEW') 
OPEN (2 1 ,FILE- S21 .RES',STATUS= r NEW') 
OPEN (22,FILE='S22.RES',STATUS= , NEW , ) 
OPEN (23, FILE- S23 RES', STATUS='NEW’) 
OPEN (24,FILE=’S24.RES , ,STATUS= , NEW’) 
OPEN (3 1 , FILE- S3 1 .RES , ,STATUS='NEW') 
OPEN (32,FILE- S32.RES',STATUS='NEW’) 
OPEN (33,FILE =I S33.RES’,STATUS= , NEW’) 
OPEN (34,FILE= , S34.RES',STATUS= r NEW’) 
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OPEN (4 1 ,FILE-C4 1 RES',STATUS= , NEW') 
OPEN (42, FILE- C42. RES', STATUS =r NEW') 
OPEN (43, FILE- C43. RES', STATUS =r NEW') 



READ DATA INTERACTIVELY 

WRITE (M001) 

READ (*,*) ICON 
IF (ICON EQ.2) GO TO 331 
WRITE (*,1009) 

READ (*,*) IVAR 
GO TO (33 1,332), IVAR 

331 WRITE (*,1002) 

READ (*,*) UMIN 
WRITE (*,1003) 

READ (*,*) UMAX 
WRITE (*,1004) 

READ (*,*) ITER 
JU=ITER 

IF (ICON EQ 2) GO TO 332 
WRITE (*,1010) 

READ (*,*) XG 
XG=XG *LENGTH 
GO TO 300 

332 WRITE (*,1005) 

READ (*,*) XGMIN 
WRITE (*,1006) 

READ (*,*) XGMAX 
WRITE (*,1004) 

READ (*,*) ITER 
JXG=ITER 

XGMIN=XGMIN * LEN GTH 
XGMAX=XGMAX*LENGTH 
IF (ICON. EQ.2) GO TO 334 
WHITE (*,1012) 

READ (*,*) U 
334 CONTINUE 
GO TO 300 
300 WRITE (*,1013) 

READ (*,*) I SET 
IF (ICON. EQ.2) GO TO 370 
GOTO (350,360,400), ISET 
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c 

C EXACT SOLUTION SET 
C 

350 DO 2 1=1, ITER 

WHITE (*,2000) I, ITER 
GO TO (351,352), IVAR 

351 U=UMIN+(UMAX-UMIN)* (I- 1 )/(ITER- 1 ) 

GO TO 355 

352 XG=XGMIN+(XGMAX-XGMIN)*(I- 1 )/(ITER- 1 ) 

355 A1=ZW*MD-MW*ZD 

A2=-CDZ* AREA*(MD-XA*ZD) 

A3=-(XG-XB)*WEIGHT*ZD 
A4= (ZG-ZB)*WEIGHT*ZD 
CALL 

EXS0LS(IVAR,LL,IFC,U,XG,A1,A2,A3,A4,CDZ,AREA,ZW,ZD,ZG,MW, 
& MD.BUO,XA) 

IF (LLGT 1)IFC=1 

2 CONTINUE 
GO TO 5000 

C 

C PITCHFORK SOLUTION SET 
C 

360 DO 3 1=1, ITER 

WHITE (*,2000) I, ITER 
GOTO (361,362), IVAR 

361 U=UMIN+(UMAX-UMIN)*(I-1)/(ITER-1) 

GO TO 365 

362 XG=XGMIN+(XGMAX-XGMIN)*(I-1)/(ITER- 1 ) 

365 A1=ZD*(ZG-ZB)* WEIGHT 

A2= 2 0*U*U*U*CDZ*AREA*(MD-XA*ZD) 
A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4= 2 0*U*U*ZD*(XG-XB)*WEIGHT 
A5=-ZD*(XG-XB)*WEIGHT 

CALL PITSLS(IVAR,LL,IFC,U,XG,A1,A2,A3,A4,A5) 

IF (LL GT 1 ) IFC= 1 

3 CONTINUE 
GO TO 5000 

C 

C SIMPLE PITCHFORK SOLUTION SET 
C 

400 DO 401 1=1, ITER 

WHITE (*,2000) I,ITER 
GO TO (402,403), IVAR 



100 



non 



402 U=UMIN+(UMAX-UMTN)*(I- 1 )/(ITER- 1 ) 

GO TO 405 

403 XG=XGMIN+(XGMAX-XGMIN)*(I- 1 )/(ITER- 1 ) 

405 Al = ZD*(ZG-ZB)* WEIGHT 

A3 =-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)* WEIGHT) 
A4= 2. 0*U*U*ZD*(XG-XB)* WEIGHT 
CALL SIMSLS(IVAR,LL,IFC,U,XG,A1,A3,A4) 

IF (LL GT.l) IFC=1 
401 CONTINUE 
GO TO 5000 

370 GO TO (380,390,410,430), ISET 
EXACT BIFURCATION SET 

380 DO 4 1=1, JXG 

WRITE (*,2000) I, JXG 

XG=XGMIN+(XGMAX-XGMIN)*(I- 1 )/( JXG- 1 ) 

IHS=0 

DO 5 J=1,JU 

U=UMIN+(UMAX-UMIN)*(J- 1 )/(JU- 1 ) 

Al= ZW*MD-MW*ZD 
A2=-CDZ*AREA*(MD-XA*ZD) 

A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
A4= (ZG*WEIGHT-ZB*BUO)*ZD 
CALL EXNUMS(NUM,A1,A2,A3,A4,U) 

IF (J.NE.l) GO TO 381 
UO=U 

NUMO=NUM 
GO TO 5 

381 UN=U 
NUMN=NUM 

NDIF=I AB S(NUMO-NUMN) 

IF (NDIF.EQ 2) GO TO 382 
UO=UN 
NUMO=NUMN 
GO TO 5 

382 Ul=UO 
U2=UN 

DO 61 JJ=1,JU 
U=U 1 +(U2-U 1 )*(JJ- 1 )/(JU- 1 ) 

Al= ZW*MD-MW*ZD 
A2=-CDZ*AREA*(MD-XA*ZD) 

A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
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A4= (ZG*WEIGHT-ZB*BUO)*ZD 
CALL EXNUMS(NUM,A1,A2,A3,A4,U) 

IF (JJNE.l) GO TO 383 
UO=U 

NUMO=NUM 
GO TO 61 

383 UN=U 
NUMN=NUM 

NDIF=I AB S(NUMO-NUMN) 

IF (NDIF EQ 2) GO TO 384 

UO=UN 

NUMO=NUMN 

61 CONTINUE 

GO TO 5 

384 Ul=UO 
U2=UN 

DO 62 JJ=1,JU 
U=U 1 +(U2-U 1 )*( JJ- 1 )/( JU- 1 ) 

Al= ZW*MD-MW*ZD 
A2=-CDZ*AREA*(MD-XA*ZD) 

A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
A4= (ZG*WEIGHT-ZB*BUO)*ZD 
CALL EXNUMS(NUM,A1,A2,A3,A4,U) 

IF (JJ NE.l) GO TO 386 
UO=U 

NUMO=NUM 
GO TO 62 
386 UN=U 

NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF (NDIF EQ 2) GO TO 385 

UO=UN 

NUMO=NUMN 

62 CONTINUE 

GOTO 5 

385 UP=(UO+UN)/2.0 
IHS=IHS+1 

UPL=SQRT(UP/(32.2*LENGTH)) 

XGL=XG/LENGTH 

IPF=10+IHS 

WRITE (IPF,100) XGL,UP 

UO=UN 

NUMO=NUMN 
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5 CONTINUE 
4 CONTINUE 
GO TO 5000 

PITCHFORK BIFURCATION SET 

390 DO 6 1=1, JXG 

WRITE (*,2000) I, JXG 

XG=XGMIN+(XGMAX-XGMIN)*(I- 1 )/(JXG- 1 ) 

IHS=0 
DO 7 J=1,JU 

U=UMIN+(UMAX-UMIN) * ( J- 1 )/( JU- 1 ) 

Al= ZD*(ZG-ZB)* WEIGHT 

A2= 2 0*U*U*U*CDZ*AREA*(MD-XA*ZD) 

A3 =-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)* WEIGHT) 
A4= 2 0*U*U*ZD*(XG-XB)* WEIGHT 
A5=-ZD*(XG-XB)* WEIGHT 
CALL PINUMS(NUM,A1,A2,A3,A4,A5,U) 

IF (J.NE.l) GO TO 391 
UO=U 

NUMO=NUM 
GO TO 7 

391 UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF (NDIF.EQ.2) GO TO 392 
UO=UN 
NUMO=NUMN 
GO TO 7 

392 Ul=UO 
U2=UN 

DO 71 JJ=1,JU 
U=U 1 +(U2-U 1 )* (JJ- 1 )/(JU- 1 ) 

Al= ZD*(ZG-ZB)* WEIGHT 

A2= 2.0*U*U*U*CDZ*AREA*(MD-XA*ZD) 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)* WEIGHT) 

A4= 2. 0*U*U*ZD*(XG-XB)* WEIGHT 

A5=-ZD*(XG-XB)* WEIGHT 

CALL PINUMS(NUM,A!,A2,A3,A4,A5,U) 

IF (JJ.NE 1) GO TO 393 
UO=U 

NUMO=NUM 
GO TO 71 
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UN=U 

NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF (NDIF.EQ.2) GO TO 394 
UO=UN 
NUMO=NUMN 

71 CONTINUE 

STOP 1 1 1 1 

394 Ul=UO 
U2=UN 

DO 72 JJ=1,JU 
U=U 1 +(U2-U 1 )*(JJ- 1 )/(JU- 1 ) 

Al= ZD*(ZG-ZB)* WEIGHT 

A2= 2.0*U*U*U*CDZ*AREA*(MD-XA*ZD) 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)* WEIGHT) 

A4= 2. 0*U*U*ZD*(XG-XB)* WEIGHT 

A5=-ZD*(XG-XB)* WEIGHT 

CALL PINUMS(NUM,A1,A2,A3,A4,A5,U) 

IF (JJ NE 1) GO TO 396 
UO=U 

NUMO=NUM 
GOTO 72 
396 UN=U 

NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF (NDIF.EQ 2) GO TO 395 

UO=UN 

NUMO=NUMN 

72 CONTINUE 

STOP 1112 

395 UP=(UO+UN)/2.0 
IHS=IHS+1 

UPL=SQRT(UP/(32.2*LENGTH)) 

XGL=XG/LENGTH 

IPF=10+IHS 

WRITE (IPF.IOO) XGL,UP 
UO=UN 
NUMO=NUMN 
7 CONTINUE 
6 CONTINUE 
GO TO 5000 
C 

C SIMPLE PITCHFORK BIFURCATION SET 
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410 DO 411 1=1, JXG 

WRITE (*,2000) I, JXG 

XG=XGMIN+(XGMAX-XGMIN)*(I- 1 )/(JXG- 1 ) 

IHS=0 

DO 412 J=1,JU 

U=UMIN+(UMAX-UMIN)*(J- 1 )/(JU- 1 ) 

Al= ZD*(ZG-ZB)* WEIGHT 

A3 =-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)* WEIGHT) 
A4= 2.0*U*U*ZD*(XG-XB)* WEIGHT 
CALL SINUMS(NUM,A1,A3,A4,U) 

IF (J.NE. 1) GO TO 413 
UO=U 

NUMO=NUM 
GO TO 412 

413 UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF (NDIF.EQ.2) GO TO 414 
UO=UN 
NUMO=NUMN 
GO TO 412 

414 Ul=UO 
U2=UN - 

DO 415 JJ=1,JU 
U=U 1 +(U2-U 1 )*(JJ- 1 )/( JU- 1 ) 

Al= ZD*(ZG-ZB)* WEIGHT 

A3 =-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)* WEIGHT) 
A4= 2. 0*U*U*ZD*(XG-XB)* WEIGHT 
CALL SINUMS(NUM,A1,A3,A4,U) 

IF (JJ.NE.l) GO TO 416 
UO=U 

NUMO=NUM 
GO TO 415 

416 UN=U 
NUMN=NUM 

NDIF=I AB S(NUMO-NUMN) 

IF (NDIF.EQ.2) GO TO 417 

UO=UN 

NUMO=NUMN 

415 - CONTINUE 

STOP 1 1 1 1 

417 Ul=UO 
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U2=UN 

DO 418 JJ=1,JU 
U=U 1 +(U2-U 1 )*(JJ- 1 )/( JU- 1 ) 

Al = ZD*(ZG-ZB)* WEIGHT 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4= 2 0*U*U*ZD*(XG-XB)* WEIGHT 
CALL SINUMS(NUM,A1,A3,A4,U) 

IF (JJ NE 1) GO TO 419 
UO=U 

NUMO=NUM 
GO TO 418 

419 UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF (NDIF.EQ 2) GO TO 420 
UO=UN 
NUMO=NUMN 
418 CONTINUE 
STOP 1112 

420 UP=(UO+UN)/2 0 
IHS=IHS+1 

UPL=SQRT(UP/(32.2*LENGTH)) 

XGL=XG/LENGTH 

IPF=10+IHS 

WRITE (IPF, 100) XGL,UP 
UO-UN 
NUMO=NUMN 
412 CONTINUE 
4 1 1 CONTINUE 
GO TO 5000 

ANALYTIC BIFURCATION SET 

0 WRITE (*,1017) 

RE.AD (*,*) ICUSP 
GOTO (432,433), ICUSP 

"EXACT" CUSP 

432 DO 431 1=1, JU 

- WHITE (*,2000) I,JU 
U=UMIN+(UMAX-UMIN)*(I- 1 )/(JU- 1 ) 
DELTA=32.2*(ZW*MD-MW*ZD)/(ZD* WEIGHT) 
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LAMBDA=U*U/(32.2*(ZG-ZB)) 

BETA2=(8.0*U*U*(1.0+DELTA*LAMBDA)**3)/(27.0*LAMBDA**2) 

IF (BETA2.LT. 0.0) GO TO 43 1 

BETA=SQRT(BETA2) 

XGB=U*U*BETA/32.2 
WHITE (10,100) U, XGB/LENGTH 
43 1 CONTINUE 
GO TO 5000 

' APPROXIMATE" CUSP 

433 DO 434 1=1, JU 

WHITE (*,2000) I,JU 

U=UMIN+(UMAX-UMIN)*(I- 1 )/(JU- 1 ) 

DELTA=32.2*(ZW*MD-MW*ZD)/(ZD*WEIGHT) 

LAMBDA=U*U/(32.2*(ZG-ZB)) 

BETA2=(8.0*U*U*(1.0+DELTA*LAMBDA)**3)/27.0 

IF (BETA2.LT. 0.0) GO TO 434 

BETA=SQRT(BETA2) 

XGB=U*U*BETA/32.2 
WHITE (10,100) U, XGB/LENGTH 

434 CONTINUE 
C 

5000 STOP 

1001 FORMAT (’ ENTER 1 : SOLUTION SETS ’,/ 

1 '2 BIFURCATION GRAPHS ') 

1009 FORMAT (' ENTER 1 : U VARIATION ',/ 

1 ' 2 : XG VARIATION ') 

1002 FORMAT (' ENTER MINIMUM VALUE OF U’) 

1003 FORMAT (' ENTER MAXIMUM VALUE OF U') 

1004 FORMAT (' ENTER NUMBER OF INCREMENTS') 

1005 FORMAT (' ENTER MINIMUM VALUE OF XG/L') 

1006 FORMAT (’ ENTER MAXIMUM VALUE OF XG/L') 

1010 FORMAT 0 ENTER VALUE OF XG/L') 

1012 FORMAT (' ENTER VALUE OF U') 

1013 FORMAT (' ENTER 1 : EXACT COMPUTATION ',/ 

1 ' 2 : PITCHFORK APPROXIMATION',/ 

2 ' 3 . SIMPLE PITCHFORK ',/ 

3 ' 4 : ANALYTIC SET ') 

1017 FORMAT (' ENTER 1 : EXACT CUSP ',/ 

1 ' 2 : APPROXIMATE CUSP ') 

2000 FORMAT (215) 

100 FORMAT (2E15.5) 
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END 



SUBROUTINE TRAP(N,A,B,OUT) 

NUMERICAL INTEGRATION ROUTINE USING THE TRAPEZOIDAL RULE 

DIMENSION A(1),B(1) 

N1=N-1 
OUT=0 0 
DO 1 I=1,N1 

OUT 1 =0. 5 *( A(I)+ A(I+ 1 ))*(B(I+ 1 )-B(I)) 

OUT =OUT+OUTl 
1 CONTINUE 
RETURN 
END 



SUBROUTINE EXSOLS(IVAR,L,IFC,U,XG,AI,A2,A3,A4,CDZ,AREA,ZW, 
&ZD.ZG,MW,MD,BUO,XA) 

IT COMPUTES AND PRINTS THE EXACT SOLUTION SET THETA 
(DEG ) VERSUS U OR XG 

RE AL MW,MD,PRNT,PI 
DIMENSION VF(5,2) 

PI=4 0*ATAN(1.0) 

IF (IVAR.EQ. 1) PRNT=U 
IF (IVAR.EQ. 2) PRNT=XG 

FIND FIRST ESTIMATE OF SOLUTIONS 



L=0 

VA=-89.5 

VA=VA*PI/180.0 

VAO=VA 

VO=THETEQ( 1 , VA, A 1 , A2, A3, A4,U) 
VAMIN=-89.5 
VAMAX=+89 5 
IVA=5000 
DO 10 I=2,IVA 
AI=I 

V A=V AMIN+C V AMAX-VAMIN)*( I- 1 )/(I VA- 1 ) 
C VA=AI-90.5 

VA=VA*PI/1 80.0 
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VAN=VA 

VN=THETEQ( 1 , VA, A 1 , A2, A3 , A4, U) 

VP=VO*VN 

IF (VP GE.0.0) GO TO 11 
L=L+1 

VF(L,1)=VA0 
VF(L,2)=VAN 
1 1 VO=VN 
VAO=VAN 
10 CONTINUE 

EXACT COMPUTATION OF SOLUTIONS VIA NEWTON’S METHOD 

E=l.E-8 
IEND=500 
DO 20 J=1,L 

X=( VF(J, 1 )+VF(J,2))/2.0 
C 

X1=X 

IXX=0 

IF (IXX EQ.O) GO TO 35 
C 

F=THETEQ(1,X,A1,A2,A3,A4,U) 

FDER=THETEQ(2,X,A1,A2,A3,A4,U) 

DO 30 K=1,IEND 
IF (FDER.EQ.0.0) STOP 1001 
DX=F/FDER 
X1=X-DX 

F=THETEQ( 1 ,X 1 , A1 , A2, A3, A4,U) 

FDER=THETEQ(2, X 1 , A 1 , A2, A3 , A4, U) 

IF (F.EQ.0.0) GO TO 35 
A=ABS(X1-X) 

IF (A-E) 35,35,40 
40 X=X1 
30 CONTINUE 
GOTO 20 

35 SOL=Xl* 180/PI 
C WRITE (*,*) XI 
JPR=10+J 
WD=U*TAN(X1) 

DLT=((CDZ*AREA*WD : "ABS(WD))-(ZW*U*WD))/(ZD*U*U) 

IF ((ABS(DLT).GT.0.4).AND.(CDZ.NE.0.0)) THEN 
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CALL 

SATUR(U,XG,CDZ,AREA,ZW,ZD,ZG,MW,MD,BUO,XA,DLT,PRNT) 

ELSELF ((ABS(DLT) GT.0.4) AND.(CDZ.EQ.O.O)) THEN 
CALL 

SATUR1 (U,XG,CDZ, ARE A,ZW,ZD,ZG,MW,MD,BUO,XA,DLT, PRNT) 

END IF 

IF (L EQ.l .AND.IFC EQ 0) WRITE (10,100) PRNT, SOL, DLT* 1 80/PI 
IF (L EQ LAND IFC EQ 1) WRITE (16,100) PRNT,SOL,DLT* 180/PI 
IF (L.GT.l) WRITE (JPRJOO) PRNT,SOL,DLT* 1 80/PI 

20 CONTINUE 
RETURN 

100 FORMAT (3E15.5) 

END 

C 

SUBROUTINE 

SATUR(U,XG,CDZ, ARE A,ZW,ZD,ZG,MW,MD,BUO,XA,DLT, PRNT) 

REAL MW, MD, PRNT 
C 

C USED TO COMPUTE SATURATION CONDITIONS FOR THE DIVE PLANE 
WITH DRAG 

C COEFFICENT INCLUDED 
C 

CD=CDZ 

IF (DLT GT 0 0) THEN 

Wl=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U*0.4)) 
& /(2*(-CD*AREA)) 

IF (W1.GT.0.0) THEN 

A=MW*U*W1-CD*XA*AREA*W1*W1-XG*BUO+MD*U*U*0 4 
B=-2.0*ZG*BUO 

C=MW*U*W1 -CD*XA* AREA*W1 *W1 +XG*BUO+MD*U*U*0.4 
X 1 =(-B+SQRT((B 1 *B 1 )-4*(A*C)))/(2.0*A) 

SOL1=2.0*ATAN(X1) 

X2=(-B-SQRT((B 1 *B 1 )-4*(A*C)))/(2.0*A) 

SOL2=2 0*ATAN(X2) 

WRITE (21,200) PRNT,SOL 1 ,SOL2 
END IF 

W2=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U*0.4)) 
& /(2* (-CD* AREA)) 

IF (W2 GT.0.0) THEN 

A=MW*U*W2-CD*XA*AREA*W2*W2-XG*BUO+MD*U*U*0 4 
B=-2.0*ZG*BUO 

C=MW*U*W2-CD*XA*AREA*W2*W2+XG*BUO+MD*U*U*0.4 
X1=(-B+SQRT((B1 *B1)-4*(A*C)))/(2.0*A) 
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SOL1=2.0*ATAN(X1) 

X2=(-B-SQRT((B 1 *B 1 )-4*(A*C)))/(2 0*A) 

SOL2=2.0*ATAN(X2) 

WRITE (22,200) PRNT,SOLl,SOL2 
END IF 

W3=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U*0 4)) 
& /(2*(CD*AREA)) 

EF (W3.LT 0.0) THEN 

A=MW*U*W3+CD*XA*AREA*W3*W3-XG*BUO+MD*U*U*0.4 
B=-2 0*ZG*BUO 

C=MW*U*W3+CD*XA*AREA*W3*W3+XG*BUO+MD*U*U*0 4 

X 1=(-B+SQRT((B 1 *B 1 )-4*(A*C)))/(2.0* A) 

SOL1=2.0*ATAN(X1) 

X2=(-B-SQRT((B 1 *B1)-4*(A*C)))/(2.0*A) 

SOL2=2.0*ATAN(X2) 

WRITE (23,200) PRNT,SOLl,SOL2 
END IF 

W4=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U*0 4)) 
& /(2*(CD*AREA)) 

IF (W4.LT.O.O) THEN 

A=MW*U*W4+CD*XA*AREA*W4*W4-XG*BUO+MD*U*U*0.4 
B=-2 0*ZG*BUO 

C=MW*U*W4+CD*XA* AREA* W4*W4+XG*BUO+MD*U*U *0.4 
X1=(-B+SQRT((B1*B1)-4*(A*C)))/(2.0*A) 

SOL1=2.0*ATAN(XI) 

X2=(-B-SQRT((B 1 *B 1 )-4*(A*C)))/(2.0* A) 

SOL2=2.0*ATAN(X2) 

WRITE (24,200) PRNT,SOLl,SOL2 
END IF 

ELSE 

Wl=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U* 

& (-0.4)))/(2*(-CD*AREA)) 

IF (W1.GT.0.0) THEN 

A=MW*U*W1-CD*XA*AREA*W1*W1-XG*BUO+MD*U*U*(-0 4) 
B=-2.0*ZG*BUO 

C=MW*U*W1-CD*XA*AREA*W1*W1+XG*BUO+MD*U*U*(-0 4) 

X I =(-B+SQRT((B 1 *B 1 )-4*( A*C)))/(2 0* A) 

SOL 1=2.0* AT AN(X1) 

X2=(-B-SQRT((B 1 *B 1 )-4*(A*C)))/(2.0* A) 

SOL2=2.0*ATAN(X2) 

WRITE (3 1 ,200) PRNT,SOL 1 ,S0L2 
END IF 

W2=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U* 
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& (-0.4))) /(2*(-CD* AREA)) 

IF (W2.LT.0.0) THEN 

A=MW*U*W2-CD*XA*AREA*W2*W2-XG*BUO+MD*U*U*(-0.4) 

B=-2.0*ZG*BUO 

C=MW*U*W2-CD*XA*AREA*W2*W2+XG*BUO+MD*U*U*(-0.4) 
X 1 =(-B+SQRT ((B 1 *B 1 )-4*(A*C)))/(2 .0* A) 

SOL 1=2.0* AT AN(X1) 

X2=(-B-SQRT((B 1 *B 1 )-4*(A*C)))/(2.0* A) 

SOL2=2.0*ATAN(X2) 

WRITE (32,200) PRNT,SOL21,SOL22 
END IF 

W3=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U* 

& (-0 4)))/(2*(CD*AREA)) 

IF (W3.LT 0 0) THEN 

A=MW*U*W3+CD*XA*AREA*W3*W3-XG*BUO+MD*U*U*(-0.4) 

B=-2.0*ZG*BUO 

C=MW*U*W3+CD*XA*AREA*W3*W3+XG*BUO+MD*U*U*(- 

0 4) 

X1=(-B+SQRT((B1*B1)-4*(A*C)))/(2.0*A) 

S0L1=2 0*ATAN(X1) 

X2=(-B-SQRT((B1*B1)-4*(A*C)))/(2.0*A) 

SOL2=2.0*ATAN(X2) 

WRITE (33,200) PRNT,S0L1,S0L2 
END IF 

W4=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD* AREA)*ZD*U*U* 

& (-0.4)))/(2*(CD*AREA)) 

IF (W4 LT 0 0) THEN 

A=MW*U*W4+CD*XA*AREA*W4*W4-XG*BUO+MD*U*U*(-0.4) 
B=-2 0*ZG*BUO 

C=MW*U*W4+CD*XA*AREA*W4*W4+XG*BUO+MD*U*U*(- 

0 4) 

X1=(-B+SQRT((B1 *B1)-4*(A*C)))/(2.0*A) 

SOL 1=2.0* AT AN(X1) 

X2=(-B-SQRT((B 1 *B 1 )-4*(A*C)))/(2 0* A) 

SOL2=2.0*ATAN(X2) 

WRITE (34,200) PRNT,SOLl,SOL2 
END IF 
END IF 

200 FORMAT (3E15. 5) 

END 

C 

SUBROUTINE 

SATUR1(U,XG,CDZ,AREA,ZW,ZD,ZG,MW,MD,BU0,XA,DLT,PRNT) 
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REAL MW,MD,PRNT 



COMPUTES DIVE PLANE SATURATION NEGLECTING DRAG 
PI=4.0*ATAN(1.0) 

A1=MD*U*U*0.4-MW*U*U*ZD*0 4-XG*ZW*BU0 
A2=-MD*U*U*0.4+MW*U*U*ZD*0.4-XG*ZW*BUO 
Bl=-2 0*ZG*ZW*BUO 

C1=-MW*U*U*ZD*0.4+XG*ZW*BUO+MD*U*U*ZW*0 4 

C2=MW*U*U*ZD*0.4+XG*ZW*BUO-MD*U*U*ZW*0.4 

XCDl=((-Bl/Al)+SQRT(((Bl/Al)**2)-4.0*(Cl/Al)))/2.0 

XCD2=((-B 1 /A 1 )-SQRT(((B 1/A 1 )* *2)-4.0*(C 1/A 1 )))/2.0 

XCD3=((-Bl/A2)+SQRT(((Bl/A2)**2)-4 0*(C2/A2)))/2 0 

XCD4=((-B 1 /A2)-SQRT(((B 1/A2)* *2)-4.0*(C2/A2)))/2.0 

XC 1 =2.0* ATAN(XCDl)* 180/PI 

XC2=2.0*ATAN(XCD2)* 1 80/PI 

XC3=2 0*ATAN(XCD3)* 1 80/PI 

XC4=2 0*ATAN(XCD4)* 1 80/PI 

WRITE (20, 1 00) PRNT,XC 1 ,DLT* 1 80/PI 

WRITE (41,1 00) PRNT,XC2,DLT* 1 80/PI 

WRITE (42,100) PRNT,XC3,DLT* 180/PI 

WRITE (43, 1 00) PRNT,XC4,DLT* 1 80/PI 

0 FORMAT (3E 15 5) 

END 



SUBROUTINE PITSLS(IVAR,LL,IFC,U,XG,A1,A2,A3,A4,A5) 

IT COMPUTES AND PRINTS THE PITCHFORK SOLUTION SET THETA 
(DEG.) VERSUS U OR XG 

DIMENSION VF(5,2) 

PI=4.0*ATAN(1 .0) 

IF (IVAR.EQ 1) PRNT=U 
IF (IVAR.EQ. 2) PRNT=XG 

FIND FIRST ESTIMATE OF SOLUTIONS 



L=0 

VA=-89.5 

VA=VA*PI/180.0 

VAO=VA 

VO=PITCEQ( 1 , VA, A1 , A2,A3, A4, A5,U) 
DO 10 1=2,180 
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AI=I 

VA=AI-90 5 
VA=VA*PI/ 180.0 
VAN=VA 

VN=PITCEQ(1,VA,A1,A2,A3,A4,A5,U) 

VP=VO*VN 

IF (VP GE O 0) GO TO 11 
L=L+1 

VF(L,l)=VAO 
VF(L,2)=VAN 
1 1 VO=VN 
VAO=VAN 
10 CONTINUE 

EXACT COMPUTATION OF SOLUTIONS VIA NEWTON'S METHOD 

E=1 E-5 
IEND=500 
DO 20 J=1,L 

X=( VF(J, 1 )+ VF(J,2))/2.0 
F=PITCEQ(1,X,A1,A2,A3,A4,A5,U) 

X1=X 

IXX-0 

IF (IXX EQ 0) GO TO 35 

FDER=PITCEQ(2,X,A1,A2,A3,A4,A5,U) 

DO 30 K=1,IEND 
IF (FDER EQ 0.0) STOP 1001 
DX=F/FDER 
X1=X-DX 

F=PITCEQ(1,X1,A1,A2,A3,A4,A5,U) 

FDER=PITCEQ(2,X1,A1,A2,A3,A4,A5,U) 

IF (F.EQ.0.0) GO TO 35 
A=ABS(X1-X) 

IF (A-E) 35,35,40 
40 X=X1 
30 CONTINUE 
GO TO 20 

35 SOL=Xl* 180.0/PI 
JPR=10+J 

IF (L.EQ.l .AND EFC.EQ 0) WRITE (10,100) PRNT.SOL 
IF (L.EQ.l.AND IFC EQ 1) WRITE (16,100) PRNT,SOL 
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WRITE ( JPR, 1 00) PRNT,SOL 



IF(LGT.l) 

20 CONTINUE 
RETURN 

100 FORMAT (2E15.5) 
END 



C 

SUBROUTINE SIMSLS(IVAR,L,IFC,U,XG,A1,A3,A4) 



IT COMPUTES AND PRINTS THE SIMPLE PITCHFORK SOLUTION SET 



(DEG.) VERSUS U OR XG 

DIMENSION VF(5,2) 

PI=4.0*ATAN( 1 .0) 

IF (IVAR.EQ. 1 ) PRNT=U 
IF (IVAR.EQ.2) PRNT=XG 

FIND FIRST ESTIMATE OF SOLUTIONS 



L=0 

VA=-89.5 
VA=VA*PI/1 80.0 
VAO=VA 

VO=SEMPEQ( 1 , VA, A 1 , A3 , A4,U) 

DO 10 1=2,180 
AI=I 

VA=AI-90.5 

VA=VA*PI/180.0 

VAN=VA 

VN=SIMPEQ(1,VA,A1,A3,A4,U) 

VP=VO*VN 

IF (VP.GE.0.0) GO TO 1 1 
L=L+1 

VF(L, l)=VAO 
VF(L,2)=VAN 
11 VO=VN 
VAO=VAN 
10 CONTINUE 

EXACT COMPUTATION OF SOLUTIONS VIA NEWTON'S METHOD 

E=l.E-5 

IEND=500 
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DO 20 J=1,L 

X=(VF(J, 1 )+VF(J,2))/2.0 
C 

X1=X 

IXX=0 

IF (IXX EQ.O) GO TO 35 
C 

F=SIMPEQ( 1 ,X, A 1 , A3, A4,U) 

FDER=SEMPEQ(2,X,A1,A3,A4,U) 

DO 30K=1,IEND 
IF (FDER.EQ 0 0) STOP 1001 
DX=F/FDER 
X1=X-DX 

F=SIMPEQ(1 ,X1,A1,A3,A4,U) 

FDER=SEMPEQ(2,X1 ,A1,A3,A4,U) 

IF (F.EQ.0.0) GO TO 35 
A=.ABS(X1-X) 

IF (A-E) 35,35,40 
40 X=X1 
30 CONTINUE 
GOTO 20 

35 SOL=Xl* 180.0/PI 
JPR= 1 0+J 

IF (L EQ 1 AND IFC.EQ 0) WRITE (10,100) PRNT,SOL 
IF (L EQ 1. AND. IFC.EQ 1) WRITE (16,100) PRNT,SOL 
IF(L.GT l) WRITE ( JPR,1 00) PRNT, SOL 

20 CONTINUE 
RETURN 

100 FORMAT (2E15.5) 

END 

C 

SUBROUTINE EXNUMS(NUM,A1,A2,A3,A4,U) 

IT COMPUTES THE NUMBER OF SOLUTIONS OF THE EXACT SOLUTION 

SET 

C 

PI=4.0*ATAN( 1 .0) 

L=0 

VA=-89.5 

VA=VA*PI/180,0 

VO-THETEQ( 1 , VA, A1 , A2, A3, A4,U) 

DO 10 1=2,180 
AI=I 
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VA=AI-90.5 

VA=AI*PI/180.0 

VN=THETEQ(1,VA,A1,A2,A3,A4,U) 

VP=VO*VN 

IF (VP.GE.0.0) GO TO 1 1 
L=L+1 
1 1 VO=VN 
10 CONTINUE 
NUM=L 
RETURN 
END 

C 

SUBROUTINE PINUMS(NUM,A1,A2,A3,A4,A5,U) 



IT COMPUTES THE NUMBER OF SOLUTIONS OF THE PITCHFORK 
OLUTION SET 



PI=4.0*ATAN(1.0) 

L=0 

VA=-89 5 
VA=VA*PI/180.0 

VO=PITCEQ( 1 , VA, A 1 , A2, A3 , A4, A5,U) 

DO 10 1=2,180 
AI=I 

VA=AI-90.5 

VA=VA*PI/180.0 

VN=PITCEQ(1,VA,A1,A2,A3,A4,A5,U) 

VP=VO*VN 

IF (VP.GE.0.0) GO TO 11 
L=L+1 
11 VO=VN 
10 CONTINUE 
NUM=L 
RETURN 
END 



SUBROUTINE SINUMS(NUM,A1,A3,A4,U) 

IT COMPUTES THE NUMBER OF SOLUTIONS OF THE SIMPLE PITCHFORK 
SOLUTION SET 

PI=4.0* ATAN( 1 .0) 

L=0 
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VA=-89.5 

VA=VA*PI/180 0 

V0=SIMPEQ(1,VA,A1,A3,A4,U) 

DO 10 1=2,180 
AI=I 

VA=AI-90.5 

VA=VA*PI/180.0 

VN=SEMPEQ(1,VA,A1,A3,A4,U) 

VP=VO*VN 

IF (VP GE O 0) GO TO 1 1 
L=L+1 
11 VO=VN 
10 CONTINUE 
NUM=L 
RETURN 
END 

C 

FUNCTION THETEQ(K,XA,A1,A2,A3,A4,U) 

C 

C IT COMPUTES THE VALUE OF THE EXACT EQUATION FOR 
C THETABAR FOR A GIVEN VALUE OF THETA 
C 

C K = 1 : COMPUTE THE VALUE OF THE FUNCTION 
C K = 2 : COMPUTE THE VALUE OF ITS DERIVATIVE 
C 

SP=SIN(XA) 

CP=COS(XA) 

AP=ABS(SP) 

TP=TAN(XA) 

ATP=ABS(TP) 

GO TO (10,20), K 

10 THETEQ=A1 *U*U*SP*CP+A2*U*U*SP*AP+A3*CP**3+A4*SP*CP**2 
C 10 THETEQ=A1 *U*TP+A2*U*U*TP*ATP+A3*CP+A4*SP 
GOTO 50 

20 THETEQ=A1*U*U*CP*CP-A1*U*U*SP*SP-3.0*A3*CP*CP*SP+A4*CP**3 
& -2.0*A4*SP*SP*CP+2.0*A2*U*U*CP*AP 

50 RETURN 
END 



FUNCTION PITCEQ(K,XA,A1,A2,A3,A4,A5,U) 

C 

C IT COMPUTES THE VALUE OF THE PITCHFORK EQUATION FOR 
C THETABAR FOR A GIVEN VALUE OF THETA 
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K = 1 : COMPUTE THE VALUE OF THE FUNCTION 
K = 2 COMPUTE THE VALUE OF ITS DERIVATIVE 

TP =TAN(XA) 

CP =COS(XA) 

CP2=CP*CP 
AP =ABS(TP) 

GO TO (10,20), K 

10PITCEQ=A1*(U*TP)**3+A2*U*U*TP*AP+A3*U*TP+A4+A5*U*U*TP*TP 
GOTO 50 

20 PITCEQ=3.0*U*U*U*TP/CP2+2 0*A2*U*U*AP/CP2+A3*U/CP2 
& +2.0*A5*U*U*TP/CP2 

50 RETURN 
END 



FUNCTION SIMPEQ(K,XA,A1,A3,A4,U) 

C 

C IT COMPUTES THE VALUE OF THE SIMPLE PITCHFORK EQUATION FOR 
C THETABAR FOR A GIVEN VALUE OF THETA 
C 

C K = 1 : COMPUTE THE VALUE OF THE FUNCTION 
C K = 2 : COMPUTE THE VALUE OF ITS DERIVATIVE 
C 

TP =TAN(XA) 

CP =COS(XA) 

CP2=CP*CP 
AP =ABS(TP) 

GO TO (10,20), K 

10 SEMPEQ=A1 *(U*TP)**3+A3 *U*TP+A4 
GOTO 50 

20 SIMPEQ=3 0*U*U*U*TP/CP2+A3 *U/CP2 
50 RETURN 
END 
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% THIS IS PROGRAM BEFSU M. IT COMPUTES THE BIASED BIFURCATION 
CUSP FOR 

% DELTA SATURATED INCLUDING DRAG TERMS. IT USES THE FZERO 
FUNCTION TO 
% SOLVE THE EQUATIONS. 



global WBlUnZw Zdlt A XA CD rho, 

% ESTABLISH GEOMETRIC PARAMETERS 
W=1 556.23 63 ;B 1 = 1 ,0001* W; 

A= 1 9. 8473 ,XA=0. 20 1 26, 

CD=0 3; 

Iy=561 32; 
g=32.2, 
m=W/g; 
rho=1.94, 

L= 13 9792, 
xb=0. 



Alphas=l, 

AJphab=0; 
zg=0. 1 , 
zb=0, 

zgb=zg-zb, 

% NON-DIMENSIONING FACTORS 

ndl=.5*rho*L A 2; 

nd2=.5*rho*L A 3; 

nd3=.5*rho*L A 4, 

nd4=.5*rho*L A 5, 

% HYDRODYNAMIC COEFFICIENTS 

Zqdnd=-6.33e-4;Zwdnd=-l 4529e-2,Zqnd=7. 545e-3;Zwnd=- 1.391 e-2, 
Zds=-5 603e-3;Zdb=0 5*(-5.603e-3);Zdltnd=(Alphas*Zds+Alphab*Zdb), 
Mqdnd=-8 . 8e-4;Mwdnd=-5.6 1 e-4;Mqnd=-3 . 702e-3 ;Mwnd= 1 .03 24e-2; 
Mds=-0 002409;Mdb=0.5*(0.002409);Mdltnd=(AIphas*Mds+AJphab*Mdb), 

Zqd=nd3 *Zqdnd;Zwd=nd2*Zwdnd,Zq=nd2*Zqnd,Zw=nd 1 *Zwnd, 

Zdlt=ndl *Zdltnd, 

Mqd=nd4*Mqdnd;Mwd=nd3*Mwdnd;Mq=nd3*Mqnd;Mw=nd2*Mwnd, 

Mdlt=nd2*Mdltnd; 
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% ESTABLISH SPEED RANGE 
U= 9:0.01:2.6; 
for n=l:length(U); 

% SOLVES FORCE EQUATION USING FZERO FUNCTION 
th(n)=fzero('satur',(-Zdlt/Zw), 0.00001); 

% COMPUTES RELATIONSHIP BETWEEN U AND Xgb WITH RESULTS FROM 
ABOVE EQUATION 

xg(n)=(Mw*U(n) A 2*tan(th(n))-zg*W*sin(th(n))-(0.5)*rho*CD*A*XA*U(n) A 2* . 

tan(th(n))*abs(tan(th(n)))+Mdlt*U(n) A 2*(0 4))/(W*cos(th(n))), 

end, 

XgbL=xg/13.9792; 

Lp=sqrt(U A 2/(3.22)), 
for n=l :length(U); 

th 1 (n)=fzero('satur 1 ’,(-Zdlt/Zw), 0.0000 1 ); 

xgl(n)=(Mw*U(n) A 2*tan(thl(n))-zg*W*sin(thl(n))-(0.5)*rho*CD*A*... 

XA*U(n) A 2*tan(thl(n))*abs(tan(thl(n)))+Mdlt*U(n) A 2*(-0.4))/(W*cos(thl(n))); 

end, 

XgbLl=xg 1/1 3.9792, 



% THIS IS PROGRAM SATUR.M. IT ESTABLISHES THE FUNCTION TO BE 
SOLVED IN 

% THE BIFSU.M PROGRAM. 



function y=satur(th); 

y=Zw*U(n) A 2*tan(th)+(W-Bl)*cos(th)+Zdlt*U(n) A 2*(0.4).. 

-(0.5)*rho*CD*A*U(n) A 2*tan(th)*abs(tan(th)); 

% THIS IS PROGRAM SATUR1.M. IT ESTABLISHES THE FUNCTION TO BE 
SOLVED IN 

% THE BIFSU.M PROGRAM. 



function y l=satur 1 (th 1 ); 

yl=Zw*U(n) A 2*tan(th 1 )+(W-B 1 )*cos(th 1 )+Zdlt*U(n) A 2*(-0.4). . . 
-(0.5)*rho*CD*A*U(n) A 2*tan(thl)*abs(tan(thl)); 
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APPENDIX B 



DTRC SUBOFF MODEL CHARACTERISTICS 
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DESIGN PARAMETERS 
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